diff --git a/BayesianTools/R/classLikelihood.R b/BayesianTools/R/classLikelihood.R index 905d9a5..887f5ac 100644 --- a/BayesianTools/R/classLikelihood.R +++ b/BayesianTools/R/classLikelihood.R @@ -8,6 +8,7 @@ #' @param sampler sampler #' @seealso \code{\link{likelihoodIidNormal}} \cr #' \code{\link{likelihoodAR1}} \cr +#' @example /inst/examples/createLikelihoodHelp.R #' @export createLikelihood <- function(likelihood, names = NULL, parallel = F, catchDuplicates=T, sampler = NULL, parallelOptions = NULL){ diff --git a/BayesianTools/R/plotSensitivityOAT.R b/BayesianTools/R/plotSensitivityOAT.R index 1c2b060..d7731e7 100644 --- a/BayesianTools/R/plotSensitivityOAT.R +++ b/BayesianTools/R/plotSensitivityOAT.R @@ -2,11 +2,18 @@ #' @author Florian Hartig #' @param bayesianSetup An object of class BayesianSetup #' @param selection indices of selected parameters -#' @param equalScale if T, y axis of all plots will have the same scale +#' @param equalScale logical, if TRUE - y axis of all plots will have the same scale #' @note This function can also be used for sensitivity analysis of an arbitrary output - just create a BayesianSetup with this output. +#' @details +#' When the scale of the parameter do not match then the plot looks like this : +#' ![](plotSensitivity-ScaleDontMatch.png "Sensitivity Plot") +#' +#' In that case, we should put 'equalScale = F'. Then the plot looks like this: +#' ![](plotSensitivity-ScaleFalse.png "Sensitivity Plot") +#' #' @example /inst/examples/plotSensitivityHelp.R #' @export -plotSensitivity <- function(bayesianSetup, selection = NULL, equalScale = T){ +plotSensitivity <- function(bayesianSetup, selection = NULL, equalScale = TRUE){ if (is.null(selection)) selection = 1:bayesianSetup$numPars n = length(selection) @@ -39,7 +46,7 @@ plotSensitivity <- function(bayesianSetup, selection = NULL, equalScale = T){ for (i in 1:n){ - if(equalScale == T) plot(resp~par, xlab = names[i], type = "l", col = "red", data = post[[i]], ylim = c(minR, maxR), ylab = "Response") + if(equalScale == TRUE) plot(resp~par, xlab = names[i], type = "l", col = "red", data = post[[i]], ylim = c(minR, maxR), ylab = "Response") else plot(resp~par, xlab = names[i], type = "l", col = "red", data = post[[i]], ylab = "Response") abline(v = refPar[i]) diff --git a/BayesianTools/inst/examples/createLikelihoodHelp.R b/BayesianTools/inst/examples/createLikelihoodHelp.R new file mode 100644 index 0000000..1ccd278 --- /dev/null +++ b/BayesianTools/inst/examples/createLikelihoodHelp.R @@ -0,0 +1,63 @@ + + +data = rnorm(20) + +# create standardized likelihood density for Gaussian likelihood function +# fitting parameters mean and standard deviation + +likelihood <- function(x) { + predicted <- rep(x[1], length(data)) + LL = likelihoodIidNormal(predicted, data, x[2]) + return(LL) +} + +likelihood(x = c(1,1)) + +# The next step is optional, you can provide the likelihood density directly +# to BayesianSetup and createLikelihood will be called automatically. +# The purpose of the createLikelihood function is to parallelize the +# density function and to add some convenience functions, +# for example to catch exceptions that are raised by the call to likelihood +likelihoodClass <- createLikelihood(likelihood) + +bayesianSetup <- createBayesianSetup(likelihood = likelihoodClass, + lower = c(-10,0.001), + upper = c(10, 5)) + +settings = list(iterations = 1000) +out <- runMCMC(bayesianSetup, sampler = "DEzs", settings) +summary(out) + + + +# simulate temporally autocorrelated data +AR1sim<-function(n, a){ + x = rep(NA, n) + x[1] = 0 + for(i in 2:n){ + x[i] = a * x[i-1] + (1-a) * rnorm(1) + } + return(x) +} + +set.seed(123) +data = AR1sim(1000, 0.5) +plot(data, type = "l") + + +# create likelihood function for AR1 likelihood function +likelihood <- function(x) { + predicted <- rep(x[1], length(data)) + LL = likelihoodAR1(predicted, data, x[2], x[3]) + return(LL) +} + +bayesianSetup <- createBayesianSetup(likelihood = likelihood, + lower = c(-10,0.001, 0), + upper = c(10, 5, 0.99)) + +settings = list(iterations = 1000) +out <- runMCMC(bayesianSetup, sampler = "DEzs", settings) +summary(out, start = 200) +plot(out, start = 200) + diff --git a/BayesianTools/inst/examples/plotSensitivityHelp.R b/BayesianTools/inst/examples/plotSensitivityHelp.R index 7afbdb5..1d2109c 100644 --- a/BayesianTools/inst/examples/plotSensitivityHelp.R +++ b/BayesianTools/inst/examples/plotSensitivityHelp.R @@ -1,5 +1,22 @@ +data = rnorm(20) + +# create standardized likelihood density for Gaussian likelihood function +# fitting parameters mean and standard deviation +likelihood <- function(x) { + predicted <- rep(x[1], length(data)) + LL = likelihoodIidNormal(predicted, data, x[2]) + return(LL) +} + +setup = createBayesianSetup(likelihood, lower = c(-10, 0.01), upper = c(10,5)) + +# showing posterior response surface for parameter 1 +# note that this is plotted at +setup$prior$best +plotSensitivity(setup, selection = 2) + +plotSensitivity(setup) +# when parameters are not on the same scale, we can use +plotSensitivity(setup, equalScale = FALSE) -ll <- testDensityBanana -bayesianSetup <- createBayesianSetup(likelihood = ll, lower = rep(-10, 2), upper = rep(10, 2)) -plotSensitivity(bayesianSetup) diff --git a/BayesianTools/man/createLikelihood.Rd b/BayesianTools/man/createLikelihood.Rd index 9b9d33b..5db5de8 100644 --- a/BayesianTools/man/createLikelihood.Rd +++ b/BayesianTools/man/createLikelihood.Rd @@ -28,6 +28,71 @@ createLikelihood( } \description{ Creates a standardized likelihood class#' +} +\examples{ + + +data = rnorm(20) + +# create standardized likelihood density for Gaussian likelihood function +# fitting parameters mean and standard deviation + +likelihood <- function(x) { + predicted <- rep(x[1], length(data)) + LL = likelihoodIidNormal(predicted, data, x[2]) + return(LL) +} + +likelihood(x = c(1,1)) + +# The next step is optional, you can provide the likelihood density directly +# to BayesianSetup and createLikelihood will be called automatically. +# The purpose of the createLikelihood function is to parallelize the +# density function and to add some convenience functions, +# for example to catch exceptions that are raised by the call to likelihood +likelihoodClass <- createLikelihood(likelihood) + +bayesianSetup <- createBayesianSetup(likelihood = likelihoodClass, + lower = c(-10,0.001), + upper = c(10, 5)) + +settings = list(iterations = 1000) +out <- runMCMC(bayesianSetup, sampler = "DEzs", settings) +summary(out) + + + +# simulate temporally autocorrelated data +AR1sim<-function(n, a){ + x = rep(NA, n) + x[1] = 0 + for(i in 2:n){ + x[i] = a * x[i-1] + (1-a) * rnorm(1) + } + return(x) +} + +set.seed(123) +data = AR1sim(1000, 0.5) +plot(data, type = "l") + + +# create likelihood function for AR1 likelihood function +likelihood <- function(x) { + predicted <- rep(x[1], length(data)) + LL = likelihoodAR1(predicted, data, x[2], x[3]) + return(LL) +} + +bayesianSetup <- createBayesianSetup(likelihood = likelihood, + lower = c(-10,0.001, 0), + upper = c(10, 5, 0.99)) + +settings = list(iterations = 1000) +out <- runMCMC(bayesianSetup, sampler = "DEzs", settings) +summary(out, start = 200) +plot(out, start = 200) + } \seealso{ \code{\link{likelihoodIidNormal}} \cr diff --git a/BayesianTools/man/figures/plotSensitivity-ScaleDontMatch.png b/BayesianTools/man/figures/plotSensitivity-ScaleDontMatch.png new file mode 100644 index 0000000..fe8e932 Binary files /dev/null and b/BayesianTools/man/figures/plotSensitivity-ScaleDontMatch.png differ diff --git a/BayesianTools/man/figures/plotSensitivity-ScaleFalse.png b/BayesianTools/man/figures/plotSensitivity-ScaleFalse.png new file mode 100644 index 0000000..4d4d9a7 Binary files /dev/null and b/BayesianTools/man/figures/plotSensitivity-ScaleFalse.png differ diff --git a/BayesianTools/man/plotSensitivity.Rd b/BayesianTools/man/plotSensitivity.Rd index 4546ed4..f320320 100644 --- a/BayesianTools/man/plotSensitivity.Rd +++ b/BayesianTools/man/plotSensitivity.Rd @@ -4,27 +4,51 @@ \alias{plotSensitivity} \title{Performs a one-factor-at-a-time sensitivity analysis for the posterior of a given bayesianSetup within the prior range.} \usage{ -plotSensitivity(bayesianSetup, selection = NULL, equalScale = T) +plotSensitivity(bayesianSetup, selection = NULL, equalScale = TRUE) } \arguments{ \item{bayesianSetup}{An object of class BayesianSetup} \item{selection}{indices of selected parameters} -\item{equalScale}{if T, y axis of all plots will have the same scale} +\item{equalScale}{logical, if TRUE - y axis of all plots will have the same scale} } \description{ Performs a one-factor-at-a-time sensitivity analysis for the posterior of a given bayesianSetup within the prior range. } +\details{ +When the scale of the parameter do not match then the plot looks like this : +\figure{plotSensitivity-ScaleDontMatch.png}{Sensitivity Plot} + +In that case, we should put 'equalScale = F'. Then the plot looks like this: +\figure{plotSensitivity-ScaleFalse.png}{Sensitivity Plot} +} \note{ This function can also be used for sensitivity analysis of an arbitrary output - just create a BayesianSetup with this output. } \examples{ +data = rnorm(20) + +# create standardized likelihood density for Gaussian likelihood function +# fitting parameters mean and standard deviation +likelihood <- function(x) { + predicted <- rep(x[1], length(data)) + LL = likelihoodIidNormal(predicted, data, x[2]) + return(LL) +} + +setup = createBayesianSetup(likelihood, lower = c(-10, 0.01), upper = c(10,5)) + +# showing posterior response surface for parameter 1 +# note that this is plotted at +setup$prior$best +plotSensitivity(setup, selection = 2) + +plotSensitivity(setup) +# when parameters are not on the same scale, we can use +plotSensitivity(setup, equalScale = FALSE) -ll <- testDensityBanana -bayesianSetup <- createBayesianSetup(likelihood = ll, lower = rep(-10, 2), upper = rep(10, 2)) -plotSensitivity(bayesianSetup) } \author{ Florian Hartig diff --git a/Development/Figures/figures.R b/Development/Figures/figures.R new file mode 100644 index 0000000..0ecad16 --- /dev/null +++ b/Development/Figures/figures.R @@ -0,0 +1,21 @@ +##### plotSensitivity ##### + +data = rnorm(20) + +# create standardized likelihood density for Gaussian likelihood function +# fitting parameters mean and standard deviation +likelihood <- function(x) { + predicted <- rep(x[1], length(data)) + LL = likelihoodIidNormal(predicted, data, x[2]) + return(LL) +} + +setup = createBayesianSetup(likelihood, lower = c(-10, 0.01), upper = c(10,5)) + +png(filename = "BayesianTools/man/figures/plotSensitivity-ScaleDontMatch.png", width = 400, height = 300) +plotSensitivity(setup) +dev.off() + +png(filename = "BayesianTools/man/figures/plotSensitivity-ScaleFalse.png", width = 400, height = 300) +plotSensitivity(setup, equalScale = F) +dev.off() \ No newline at end of file