diff --git a/.gitignore b/.gitignore index bd48ad3..ff72220 100644 --- a/.gitignore +++ b/.gitignore @@ -7,6 +7,7 @@ soobench/ *.html !docs/*.html +!docs/*/*.html *.gz *.dll *.DS_Store diff --git a/docs/articles/BayesianTools.html b/docs/articles/BayesianTools.html new file mode 100644 index 0000000..6330a61 --- /dev/null +++ b/docs/articles/BayesianTools.html @@ -0,0 +1,809 @@ + + +
+ + + + +Abstract
+The BayesianTools (BT) package supports model analysis (including sensitivity analysis and uncertainty analysis), Bayesian model calibration, as well as model selection and multi-model inference techniques for system models.
+The purpose of this first section is to give you a quick overview of the most important functions of the BayesianTools (BT) package. For a more detailed description, see the later sections
+If you haven’t installed the package yet, either run
+install.packages("BayesianTools")
Or follow the instructions on https://github.com/florianhartig/BayesianTools to install a development or an older version.
+Loading and citation
+library(BayesianTools)
+citation("BayesianTools")
##
+## To cite package 'BayesianTools' in publications use:
+##
+## Florian Hartig, Francesco Minunno and Stefan Paul (2017).
+## BayesianTools: General-Purpose MCMC and SMC Samplers and Tools
+## for Bayesian Statistics. R package version 0.1.4.
+## https://github.com/florianhartig/BayesianTools
+##
+## A BibTeX entry for LaTeX users is
+##
+## @Manual{,
+## title = {BayesianTools: General-Purpose MCMC and SMC Samplers and Tools for Bayesian Statistics},
+## author = {Florian Hartig and Francesco Minunno and Stefan { Paul}},
+## year = {2017},
+## note = {R package version 0.1.4},
+## url = {https://github.com/florianhartig/BayesianTools},
+## }
+Note: BayesianTools calls a number of secondary packages. Particular important is coda, which is used on a number of plots and summary statistics. If you make heavy use of the summary statistics and diagnostics plots, it would be nice to cite coda as well!
+Pro-tip: if you are running a stochastic algorithms such as an MCMC, you should always set or record your random seed to make your results reproducible (otherwise, results will change slightly every time you run the code)
+set.seed(123)
In a real application, to ensure reproducibility, it would also be useful to record the session
+sessionInfo()
which lists the version number of R and all loaded packages.
+The central object in the BT package is the BayesianSetup. This class contains the information about the model to be fit (likelihood), and the priors for the model parameters.
+A BayesianSetup is created by the createBayesianSetup function. The function expects a log-likelihood and (optional) a log-prior. It then automatically creates the posterior and various convenience functions for the samplers.
+Advantages of the BayesianSetup include 1) support for automatic parallelization, 2) functions are wrapped in try-catch statements to avoid crashes during long MCMC evaluations, 3) and the posterior checks if the parameter is outside the prior first, in which case the likelihood is not evaluated (makes the algorithms faster for slow likelihoods).
+If no prior information is provided, an unbounded flat prior is created. If no explicit prior, but lower and upper values are provided, a standard uniform prior with the respective bounds is created, including the option to sample from this prior, which is useful for SMC and also for getting starting values. This option is used in the following example, which creates a multivariate normal likelihood density and a uniform prior for 3 parameters.
+ll <- generateTestDensityMultiNormal(sigma = "no correlation")
+bayesianSetup = createBayesianSetup(likelihood = ll, lower = rep(-10, 3), upper = rep(10, 3))
See later more detailed description about the BayesianSetup.
+Hint: for an example how to run this steps for dynamic ecological model, see ?VSEM
+Once you have your setup, you may want to run a calibration. The runMCMC function is the main wrapper for all other implemented MCMC/SMC functions. It always takes the following arguments
+As an example, choosing the sampler name “Metropolis” calls a versatile Metropolis-type MCMC with options for covariance adaptation, delayed rejection, tempering and Metropolis-within-Gibbs sampling. For details, see the the later reference on MCMC samplers. This is how we would call this sampler with default settings
+iter = 10000
+settings = list(iterations = iter, message = FALSE)
+out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings)
All samplers can be plotted and summarized via the console with the standard print, and summary commands
+print(out)
## [1] "mcmcSampler - you can use the following methods to summarize, plot or reduce this class:"
+## [1] getSample plot print summary
+## see '?methods' for accessing help and source code
+summary(out)
## # # # # # # # # # # # # # # # # # # # # # # # # #
+## ## MCMC chain summary ##
+## # # # # # # # # # # # # # # # # # # # # # # # # #
+##
+## # MCMC sampler: Metropolis
+## # Nr. Chains: 1
+## # Iterations per chain: 10000
+## # Rejection rate: 0.686
+## # Effective sample size: 932
+## # Runtime: 3.06 sec.
+##
+## # Parameters
+## MAP 2.5% median 97.5%
+## par 1 0.001 -1.981 0.022 2.050
+## par 2 0.000 -2.095 -0.060 1.968
+## par 3 -0.001 -2.085 -0.048 1.967
+##
+## ## DIC: 11.72
+## ## Convergence
+## Gelman Rubin multivariate psrf: Only one chain; convergence cannot be determined!
+##
+## ## Correlations
+## par 1 par 2 par 3
+## par 1 1.000 0.004 -0.005
+## par 2 0.004 1.000 0.017
+## par 3 -0.005 0.017 1.000
+and plottted with several plot functions. The marginalPlot can either be plotted as histograms with density overlay, which is also the default, or as a violin plot (see “?marginalPlot”).
+plot(out) # plot internally calls tracePlot(out)
correlationPlot(out)
marginalPlot(out)
Other Functions that can be applied to all samplers include model selection scores such as the DIC and the marginal Likelihood (for the calculation of the Bayes factor, see later section for more details), and the Maximum Aposteriori Value (MAP). For the marginal likelihood calculation it is possible to chose from a set of methods (see “?marginalLikelihood”).
+marginalLikelihood(out)
## $ln.ML
+## [1] -8.969702
+##
+## $ln.lik.star
+## [1] -2.756817
+##
+## $ln.pi.star
+## [1] -8.987197
+##
+## $ln.pi.hat
+## [1] -2.774312
+##
+## $method
+## [1] "Chib"
+DIC(out)
## $DIC
+## [1] 11.72016
+##
+## $IC
+## [1] 14.81802
+##
+## $pD
+## [1] 3.097861
+##
+## $pV
+## [1] 3.411328
+##
+## $Dbar
+## [1] 8.622295
+##
+## $Dhat
+## [1] 5.524434
+MAP(out)
## $parametersMAP
+## par 1 par 2 par 3
+## 1.000968e-03 6.672794e-05 -9.910697e-04
+##
+## $valuesMAP
+## Lposterior Llikelihood Lprior
+## -11.744013 -2.756817 -8.987197
+You can extract (a part of) the sampled parameter values by
+getSample(out, start = 100, end = NULL, thin = 5, whichParameters = 1:2)
For all samplers, you can conveniently perform multiple runs via the nrChains argument
+iter = 1000
+settings = list(iterations = iter, nrChains = 3, message = FALSE)
+out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings)
The result is an object of mcmcSamplerList, which should allow to do everything one can do with an mcmcSampler object (with slightly different output sometimes).
+print(out)
## [1] "mcmcSamplerList - you can use the following methods to summarize, plot or reduce this class:"
+## [1] getSample plot print summary
+## see '?methods' for accessing help and source code
+summary(out)
## # # # # # # # # # # # # # # # # # # # # # # # # #
+## ## MCMC chain summary ##
+## # # # # # # # # # # # # # # # # # # # # # # # # #
+##
+## # MCMC sampler: Metropolis
+## # Nr. Chains: 3
+## # Iterations per chain: 1000
+## # Rejection rate: 0.688
+## # Effective sample size: 301
+## # Runtime: 1.11 sec.
+##
+## # Parameters
+## psf MAP 2.5% median 97.5%
+## par 1 1.004 -0.001 -1.937 -0.093 2.010
+## par 2 1.004 0.000 -1.997 0.067 1.807
+## par 3 1.027 0.000 -1.787 0.127 1.994
+##
+## ## DIC: 11.212
+## ## Convergence
+## Gelman Rubin multivariate psrf: 1.015
+##
+## ## Correlations
+## par 1 par 2 par 3
+## par 1 1.000 -0.062 -0.049
+## par 2 -0.062 1.000 0.065
+## par 3 -0.049 0.065 1.000
+For example, in the plot you now see 3 chains.
+plot(out)
There are a few additional functions that may only be available for lists, for example convergence checks
+#getSample(out, coda = F)
+gelmanDiagnostics(out, plot = T)
## Potential scale reduction factors:
+##
+## Point est. Upper C.I.
+## par 1 1.00 1.02
+## par 2 1.00 1.01
+## par 3 1.03 1.06
+##
+## Multivariate psrf
+##
+## 1.02
+The likelihood should be provided as a log density function.
+ll = logDensity(x)
See options for parallelization below. We will use a simple 3-d multivariate normal density for this demonstration.
+ll <- generateTestDensityMultiNormal(sigma = "no correlation")
+bayesianSetup = createBayesianSetup(likelihood = ll, lower = rep(-10, 3), upper = rep(10, 3))
Likelihoods are often costly to compute. If that is the case for you, you should think about parallelization possibilities. The ‘createBayesianSetup’ function has the input variable ‘parallel’, with the following options
+Algorithms in the BayesianTools package can make use of parallel computing if this option is specified in the BayesianSetup. Note that currently, parallelization is used by the following algorithms: SMC, DEzs and DREAMzs sampler. It can also be used through the BayesianSetup with the functions of the sensitivity package.
+Here some more details on the parallelization
+The in-build parallelization is the easiest way to make use of parallel computing. In the “parallel” argument you can choose the number of cores used for parallelization. Alternatively for TRUE or “auto” all available cores except for one will be used. Now the proposals are evaluated in parallel. Technically, the in-build parallelization uses an R cluster to evaluate the posterior density function. The input for the parallel function is a matrix, where each column represents a parameter and each row a proposal. In this way, the proposals can be evaluated in parallel. For sampler, where only one proposal is evaluated at a time (namely the Metropolis based algorithms as well as DE/DREAM without the zs extension), no parallelization can be used.
+The second option is to use an external parallelization. Here, a parallelization is attempted in the user defined likelihood function. To make use of external parallelization, the likelihood function needs to take a matrix of proposals and return a vector of likelihood values. In the proposal matrix each row represents one proposal, each column a parameter. Further, you need to specify the “external” parallelization in the “parallel” argument. In simplified terms the use of external parallelization uses the following steps:
+## Definition of likelihood function
+likelihood <- function(matrix){
+ # Calculate likelihood in parallel
+ # Return vector of likelihood valus
+}
+
+## Create Bayesian Setup
+BS <- createBayesianSetup(likelihood, parallel = "external" ...)
+
+## Run MCMC
+runMCMC(BS, sampler = "SMC", ...)
** Remark: even though parallelization can significantly reduce the computation time, it is not always useful because of the so-called communication overhead (computational time for distributing and retrieving infos from the parallel cores). For models with low computational cost, this procedure can take more time than the actual evaluation of the likelihood. If in doubt, make a small comparison of the runtime before starting your large sampling. **
+The prior in the BayesianSetup consists of four parts
+These information can passed by first creating an a extra object, via createPrior, or through the the createBayesianSetup function.
+You have 5 options to create a prior
+If creating a user-defined prior, the following information can/should be provided to createPrior:
+The following example from the help shows how this works
+# Create a BayesianSetup
+ll <- generateTestDensityMultiNormal(sigma = "no correlation")
+bayesianSetup = createBayesianSetup(likelihood = ll, lower = rep(-10, 3), upper = rep(10, 3))
+
+settings = list(iterations = 2500, message = FALSE)
+out <- runMCMC(bayesianSetup = bayesianSetup, settings = settings)
+
+
+newPrior = createPriorDensity(out, method = "multivariate", eps = 1e-10, lower = rep(-10, 3), upper = rep(10, 3), best = NULL)
+bayesianSetup <- createBayesianSetup(likelihood = ll, prior = newPrior)
+
+settings = list(iterations = 1000, message = FALSE)
+out <- runMCMC(bayesianSetup = bayesianSetup, settings = settings)
The runMCMC function is the central function for starting MCMC algorithms in the BayesianTools package. It requires a bayesianSetup, a choice of sampler (standard is DEzs), and optionally changes to the standard settings of the chosen sampler.
+runMCMC(bayesianSetup, sampler = “DEzs”, settings = NULL)
+One optional argument that you can always use is nrChains - the default is 1. If you choose more, the runMCMC will perform several runs.
+ll <- generateTestDensityMultiNormal(sigma = "no correlation")
+
+bayesianSetup = createBayesianSetup(likelihood = ll, lower = rep(-10, 3), upper = rep(10, 3))
+
+settings = list(iterations = 10000, nrChains= 3, message = FALSE)
+out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings)
+
+plot(out)
marginalPlot(out)
correlationPlot(out)
gelmanDiagnostics(out, plot=T)
## Potential scale reduction factors:
+##
+## Point est. Upper C.I.
+## par 1 1 1.00
+## par 2 1 1.00
+## par 3 1 1.01
+##
+## Multivariate psrf
+##
+## 1
+# option to restart the sampler
+
+settings = list(iterations = 1000, nrChains= 1, message = FALSE)
+out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings)
+
+out2 <- runMCMC(bayesianSetup = out)
+
+out3 <- runMCMC(bayesianSetup = out2)
+
+#plot(out)
+#plot(out3)
+
+# create new prior from posterior sample
+
+newPriorFromPosterior <- createPriorDensity(out2)
For convenience we define a number of iterations
+iter = 10000
The BayesianTools package is able to run a large number of Metropolis-Hastings (MH) based algorithms All of these samplers can be accessed by the “Metropolis” sampler in the runMCMC function by specifying the sampler’s settings.
+The following code gives an overview about the default settings of the MH sampler.
+applySettingsDefault(sampler = "Metropolis")
## $sampler
+## [1] "Metropolis"
+##
+## $startValue
+## NULL
+##
+## $iterations
+## [1] 10000
+##
+## $optimize
+## [1] TRUE
+##
+## $proposalGenerator
+## NULL
+##
+## $consoleUpdates
+## [1] 100
+##
+## $burnin
+## [1] 0
+##
+## $thin
+## [1] 1
+##
+## $parallel
+## NULL
+##
+## $adapt
+## [1] TRUE
+##
+## $adaptationInterval
+## [1] 500
+##
+## $adaptationNotBefore
+## [1] 3000
+##
+## $DRlevels
+## [1] 1
+##
+## $proposalScaling
+## NULL
+##
+## $adaptationDepth
+## NULL
+##
+## $temperingFunction
+## NULL
+##
+## $gibbsProbabilities
+## NULL
+##
+## $currentChain
+## [1] 1
+##
+## $message
+## [1] TRUE
+##
+## $nrChains
+## [1] 1
+##
+## $runtime
+## [1] 0
+##
+## $sessionInfo
+## R version 3.4.3 (2017-11-30)
+## Platform: x86_64-w64-mingw32/x64 (64-bit)
+## Running under: Windows 8.1 x64 (build 9600)
+##
+## Matrix products: default
+##
+## locale:
+## [1] LC_COLLATE=German_Germany.1252 LC_CTYPE=German_Germany.1252
+## [3] LC_MONETARY=German_Germany.1252 LC_NUMERIC=C
+## [5] LC_TIME=German_Germany.1252
+##
+## attached base packages:
+## [1] stats graphics grDevices utils datasets methods base
+##
+## other attached packages:
+## [1] BayesianTools_0.1.4
+##
+## loaded via a namespace (and not attached):
+## [1] Rcpp_0.12.14 DHARMa_0.1.5 knitr_1.18
+## [4] magrittr_1.5 MASS_7.3-47 tmvtnorm_1.4-10
+## [7] lattice_0.20-35 bridgesampling_0.4-0 foreach_1.4.4
+## [10] stringr_1.2.0 tools_3.4.3 grid_3.4.3
+## [13] coda_0.19-1 htmltools_0.3.6 IDPmisc_1.1.17
+## [16] iterators_1.0.9 yaml_2.1.16 rprojroot_1.3-2
+## [19] digest_0.6.13 gmm_1.6-1 numDeriv_2016.8-1
+## [22] Brobdingnag_1.2-4 Matrix_1.2-12 codetools_0.2-15
+## [25] evaluate_0.10.1 rmarkdown_1.8 sandwich_2.4-0
+## [28] stringi_1.1.6 compiler_3.4.3 backports_1.1.2
+## [31] stats4_3.4.3 mvtnorm_1.0-6 zoo_1.8-1
+The following examples show how the different settings can be used. As you will see different options can be activated singly or in combination.
+The following settings will run the standard Metropolis Hastings MCMC.
+Refernences: Hastings, W. K. (1970). Monte carlo sampling methods using markov chains and their applications. Biometrika 57 (1), 97-109.
+Metropolis, N., A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller (1953). Equation of state calculations by fast computing machines. The journal of chemical physics 21 (6), 1087 - 1092.
+settings <- list(iterations = iter, adapt = F, DRlevels = 1, gibbsProbabilities = NULL, temperingFunction = NULL, optimize = F, message = FALSE)
+out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings)
+plot(out)
This sampler uses an optimization step prior to the sampling process. The optimization aims at improving the starting values and the covariance of the proposal distribution.
+settings <- list(iterations = iter, adapt = F, DRlevels = 1, gibbsProbabilities = NULL, temperingFunction = NULL, optimize = T, message = FALSE)
+out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings)
+plot(out)
In the adaptive Metropolis sampler (AM) the information already acquired in the sampling process is used to improve (or adapt) the proposal function. In the BayesianTools package the history of the chain is used to adapt the covariance of the propoasal distribution.
+References: Haario, H., E. Saksman, and J. Tamminen (2001). An adaptive metropolis algorithm. Bernoulli , 223-242.
+settings <- list(iterations = iter, adapt = T, DRlevels = 1, gibbsProbabilities = NULL, temperingFunction = NULL, optimize = T, message = FALSE)
+out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings)
+plot(out)
Even though rejection is an essential step of a MCMC algorithm it can also mean that the proposal distribution is (locally) badly tuned to the target distribution. In a delayed rejection (DR) sampler a second (or third, etc.) proposal is made before rejection. This proposal is usually drawn from a different distribution, allowing for a greater flexibility of the sampler. In the BayesianTools package the number of delayed rejection steps as well as the scaling of the proposals can be determined. ** Note that the current version only supports two delayed rejection steps. **
+References: Green, Peter J., and Antonietta Mira. “Delayed rejection in reversible jump Metropolis-Hastings.” Biometrika (2001): 1035-1053.
+settings <- list(iterations = iter, adapt = F, DRlevels = 2, gibbsProbabilities = NULL, temperingFunction = NULL, optimize = T, message = FALSE)
+out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings)
+plot(out)
The delayed rejection adaptive Metropolis (DRAM) sampler is merely a combination of the two previous sampler (DR and AM).
+References: Haario, Heikki, et al. “DRAM: efficient adaptive MCMC.” Statistics and Computing 16.4 (2006): 339-354.
+settings <- list(iterations = iter, adapt = T, DRlevels = 2, gibbsProbabilities = NULL, temperingFunction = NULL, optimize = T, message = FALSE)
+out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings)
+plot(out)
To reduce the dimensions of the target function a Metropolis-within-Gibbs sampler can be run with the BayesianTools package. This means in each iteration only a subset of the parameter vector is updated. In the example below at most two (of the three) parameters are updated each step, and it is double as likely to vary one than varying two.
+** Note that currently adaptive cannot be mixed with Gibbs updating! **
+settings <- list(iterations = iter, adapt = T, DRlevels = 1, gibbsProbabilities = c(1,0.5,0), temperingFunction = NULL, optimize = T, message = FALSE)
+out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings)
+plot(out)
Simulated tempering is closely related to simulated annealing (e.g. Bélisle, 1992) in optimization algorithms. The idea of tempering is to increase the acceptance rate during burn-in. This should result in a faster initial scanning of the target function. To include this a tempering function needs to be supplied by the user. The function describes how the acceptance rate is influenced during burn-in. In the example below an exponential decline approaching 1 (= no influece on the acceptance rate)is used.
+References: Bélisle, C. J. (1992). Convergence theorems for a class of simulated annealing algorithms on rd. Journal of Applied Probability, 885–895.
+C. J. Geyer (2011) Importance sampling, simulated tempering, and umbrella sampling, in the Handbook of Markov Chain Monte Carlo, S. P. Brooks, et al (eds), Chapman & Hall/CRC.
+temperingFunction <- function(x) 5 * exp(-0.01*x) + 1
+settings <- list(iterations = iter, adapt = F, DRlevels = 1, gibbsProbabilities = c(1,1,0), temperingFunction = temperingFunction, optimize = T, message = FALSE)
+out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings)
+plot(out)
The BT package implements two versions of the differential evolution MCMC. In doubt, you should use the DEzs option.
+The first is the normal DE MCMC, corresponding to Ter Braak, Cajo JF. “A Markov Chain Monte Carlo version of the genetic algorithm Differential Evolution: easy Bayesian computing for real parameter spaces.” Statistics and Computing 16.3 (2006): 239-249. In this sampler multiple chains are run in parallel (but not in the sense of parallel computing). The main diference to the Metrpolis based algorithms is the creation of the propsal. Generally all samplers use the current positin of the chain and add a step in the parameter space to generate a new proposal. Whereas in the Metropolis based sampler this step is usually drawn from a multivariate normal distribution (yet every distribution is possible), the DE sampler uses the current position of two other chains to generate the step for each chain. For sucessful sampling at least 2*d chains, with d being the number of parameters, need to be run in parallel.
+settings <- list(iterations = iter, message = FALSE)
+out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "DE", settings = settings)
+plot(out)
The second is the Differential Evolution MCMC with snooker update and sampling from past states, corresponding to ter Braak, Cajo JF, and Jasper A. Vrugt. “Differential evolution Markov chain with snooker updater and fewer chains.” Statistics and Computing 18.4 (2008): 435-446. This extension covers two differences to the normal DE MCMC. First a snooker update is used based on a user defined probability. Second also past states of other chains are respected in the creation of the proposal. These extensions allow for fewer chains (i.e. 3 chains are usually enough for up to 200 parameters) and parallel computing as the current position of each chain is only dependent on the past states of the other chains.
+settings <- list(iterations = iter, message = FALSE)
+out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "DEzs", settings = settings)
+plot(out)
Also for the DREAM sampler, there are two versions included. First of all, the standard DREAM sampler, see Vrugt, Jasper A., et al. “Accelerating Markov chain Monte Carlo simulation by differential evolution with self-adaptive randomized subspace sampling.” International Journal of Nonlinear Sciences and Numerical Simulation 10.3 (2009): 273-290.
+This sampler is largely build on the DE sampler with some significant differences: 1) More than two chains can be used to generate a proposal. 2) A randomized subspace sampling can be used to enhance the efficiency for high dimensional posteriors. Each dimension is updated with a crossover probalitity CR. To speed up the exploration of the posterior DREAM adapts the distribution of CR values during burn-in to favor large jumps over small ones. 3) Outlier chains can be removed during burn-in.
+settings <- list(iterations = iter, message = FALSE)
+out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "DREAM", settings = settings)
+plot(out)
The second implementation uses the same extension as the DEzs sampler. Namely sampling from past states and a snooker update. Also here this extension allows for the use of fewer chains and parallel computing.
+Again, in doubt you should prefer “DREAMzs”.
+settings <- list(iterations = iter, message = FALSE)
+out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "DREAMzs", settings = settings)
+#plot(out)
The T-walk is a MCMC algorithm developed by Christen, J. Andrés, and Colin Fox. “A general purpose sampling algorithm for continuous distributions (the t-walk).” Bayesian Analysis 5.2 (2010): 263-281. In the sampler two independent points are used to explore the posterior space. Based on probabilities four different moves are used to generate proposals for the two points. As for the DE sampler this procedure requires no tuning of the proposal distribution for efficient sampling in complex posterior distributions.
+settings = list(iterations = iter, message = FALSE)
+
+out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Twalk", settings = settings)
All MCMCs should be checked for convergence. We recommend the standard procedure of Gelmal-Rubin. This procedure requires running several MCMCs (we recommend 3). This can be achieved either directly in the runMCMC (nrChains = 3), or, for runtime reasons, by combining the results of three independent runMCMC evaluations with nrChains = 1.
+settings <- list(iterations = iter, nrChains = 3, message = FALSE)
+out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings)
+plot(out)
#chain = getSample(out, coda = T)
+gelmanDiagnostics(out, plot = F)
## Potential scale reduction factors:
+##
+## Point est. Upper C.I.
+## par 1 1 1.00
+## par 2 1 1.01
+## par 3 1 1.01
+##
+## Multivariate psrf
+##
+## 1
+MCMCs sample the posterior space by creating a chain in parameter space. While this allows “learning” from past steps, it does not permit the parallel execution of a large number of posterior values at the same time.
+An alternative to MCMCs are particle filters, aka Sequential Monte-Carlo (SMC) algorithms. See Hartig, F.; Calabrese, J. M.; Reineking, B.; Wiegand, T. & Huth, A. Statistical inference for stochastic simulation models - theory and application Ecol. Lett., 2011, 14, 816-827
+The easiest option is to simply sample a large number of parameters and accept them according to their posterior value. This option can be emulated with the implemented SMC, setting iterations to 1.
+settings <- list(initialParticles = iter, iterations= 1)
+out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "SMC", settings = settings)
+plot(out)
The more sophisticated option is using the implemented SMC, which is basically a particle filter that applies several filter steps.
+settings <- list(initialParticles = iter, iterations= 10)
+out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "SMC", settings = settings)
+plot(out)
Note that the use of a number for initialParticles requires that the bayesianSetup includes the possibility to sample from the prior.
+There are a number of Bayesian model selection and model comparison methods. The BT implements three of the most common of them, the DIC, the WAIC, and the Bayes factor.
+On the Bayes factor, see Kass, R. E. & Raftery, A. E. Bayes Factors J. Am. Stat. Assoc., Amer Statist Assn, 1995, 90, 773-795
An overview on DIC and WAIC is given in Gelman, A.; Hwang, J. & Vehtari, A. (2014) Understanding predictive information criteria for Bayesian models. Statistics and Computing, 24, 997-1016-. On DIC, see also the original reference by Spiegelhalter, D. J.; Best, N. G.; Carlin, B. P. & van der Linde, A. (2002) Bayesian measures of model complexity and fit. J. Roy. Stat. Soc. B, 64, 583-639.
The Bayes factor relies on the calculation of marginal likelihoods, which is numerically not without problems. The BT package currently implements three methods
+The recommended way is the method “Chib” (Chib and Jeliazkov, 2001). which is based on MCMC samples, but performs additional calculations. Despite being the current recommendation, note there are some numeric issues with this algorithm that may limit reliability for larger dimensions.
The harmonic mean approximation, is implemented only for comparison. Note that the method is numerically unrealiable and usually should not be used.
The third method is simply sampling from the prior. While in principle unbiased, it will only converge for a large number of samples, and is therefore numerically inefficient.
Data linear Regression with quadratic and linear effect
+sampleSize = 30
+x <- (-(sampleSize-1)/2):((sampleSize-1)/2)
+y <- 1 * x + 1*x^2 + rnorm(n=sampleSize,mean=0,sd=10)
+plot(x,y, main="Test Data")
Likelihoods for both
+likelihood1 <- function(param){
+ pred = param[1] + param[2]*x + param[3] * x^2
+ singlelikelihoods = dnorm(y, mean = pred, sd = 1/(param[4]^2), log = T)
+ return(sum(singlelikelihoods))
+}
+
+likelihood2 <- function(param){
+ pred = param[1] + param[2]*x
+ singlelikelihoods = dnorm(y, mean = pred, sd = 1/(param[3]^2), log = T)
+ return(sum(singlelikelihoods))
+}
Posterior definitions
+setUp1 <- createBayesianSetup(likelihood1, lower = c(-5,-5,-5,0.01), upper = c(5,5,5,30))
+
+setUp2 <- createBayesianSetup(likelihood2, lower = c(-5,-5,0.01), upper = c(5,5,30))
MCMC and marginal likelihood calculation
+settings = list(iterations = 15000, message = FALSE)
+out1 <- runMCMC(bayesianSetup = setUp1, sampler = "Metropolis", settings = settings)
+#tracePlot(out1, start = 5000)
+M1 = marginalLikelihood(out1)
+M1
+
+settings = list(iterations = 15000, message = FALSE)
+out2 <- runMCMC(bayesianSetup = setUp2, sampler = "Metropolis", settings = settings)
+#tracePlot(out2, start = 5000)
+M2 = marginalLikelihood(out2)
+M2
Bayes factor (need to reverse the log)
+exp(M1$marginalLikelihod - M2$marginalLikelihod)
## numeric(0)
+BF > 1 means the evidence is in favor of M1. See Kass, R. E. & Raftery, A. E. (1995) Bayes Factors. J. Am. Stat. Assoc., Amer Statist Assn, 90, 773-795.
+Note that we would have to multiply still with the model priors to arrive at Bayesian model weights.
+The Adaptive Metropolis Algorithm (Haario et al. 2001)
+ + +AM(startValue = NULL, iterations = 10000, nBI = 0, parmin = NULL, + parmax = NULL, FUN, f = 1, eps = 0)+ +
startValue | +vector with the start values for the algorithm. Can be NULL if FUN is of class BayesianSetup. In this case startValues are sampled from the prior. |
+
---|---|
iterations | +iterations to run |
+
nBI | +number of burnin |
+
parmin | +minimum values for the parameter vector or NULL if FUN is of class BayesianSetup |
+
parmax | +maximum values for the parameter vector or NULL if FUN is of class BayesianSetup |
+
FUN | +function to be sampled from or object of class bayesianSetup |
+
f | +scaling factor |
+
eps | +small number to avoid singularity |
+
Haario, Heikki, Eero Saksman, and Johanna Tamminen. "An adaptive Metropolis algorithm." Bernoulli (2001): 223-242.
+ + +A package with general-purpose MCMC and SMC samplers, as well as plots and diagnostic functions for Bayesian statistics
+ + + +A package with general-purpose MCMC and SMC samplers, as well as plots and diagnostic functions for Bayesian statistics, particularly for process-based models.
+The package contains 2 central functions, createBayesianSetup
, which creates a standardized Bayesian setup with likelihood and priors, and runMCMC
, which allows to run various MCMC and SMC samplers.
The package can of course also be used for general (non-Bayesian) target functions.
+To use the package, a first step to use createBayesianSetup
to create a BayesianSetup, which usually contains prior and likelihood densities, or in general a target function.
Those can be sampled with runMCMC
, which can call a number of general purpose Metropolis sampler, including the Metropolis
that allows to specify various popular Metropolis variants such as adaptive and/or delayed rejection Metropolis; two variants of differential evolution MCMC DE
, DEzs
, two variants of DREAM DREAM
and DREAMzs
, the Twalk
MCMC, and a Sequential Monte Carlo sampler smcSampler
.
The output of runMCMC is of class mcmcSampler / smcSampler if one run is performed, or mcmcSamplerList / smcSamplerList if several sampler are run. Various functions are available for plotting, model comparison (DIC, marginal likelihood), or to use the output as a new prior.
+For details on how to use the packgage, run vignette("BayesianTools", package="BayesianTools").
+To get the suggested citation, run citation("BayesianTools")
+To report bugs or ask for help, post a reproducible example via the BayesianTools issue tracker on GitHub.
+Acknowledgements: The creation of this package was facilicated through meetings of Cost Action FP 1304 Profound.
+ + +Differential-Evolution MCMC
+ + +DE(bayesianSetup, settings = list(startValue = NULL, iterations = 10000, f = + -2.38, burnin = 0, thin = 1, eps = 0, consoleUpdates = 100, blockUpdate = + list("none", k = NULL, h = NULL, pSel = NULL, pGroup = NULL, groupStart = + 1000, groupIntervall = 1000), currentChain = 1, message = TRUE))+ +
bayesianSetup | +a BayesianSetup with the posterior density function to be sampled from |
+
---|---|
settings | +list with parameter settings |
+
startValue | +(optional) eiter a matrix with start population, a number to define the number of chains that are run or a function that samples a starting population. |
+
iterations | +number of function evaluations. |
+
burnin | +number of iterations treated as burn-in. These iterations are not recorded in the chain. |
+
thin | +thinning parameter. Determines the interval in which values are recorded. |
+
f | +scaling factor gamma |
+
eps | +small number to avoid singularity |
+
blockUpdate | +list determining whether parameters should be updated in blocks. For possible settings see Details. |
+
message | +logical determines whether the sampler's progress should be printed |
+
For blockUpdate the first element in the list determines the type of blocking. +Possible choices are
"none" (default), no blocking of parameters
"correlation" blocking based on correlation of parameters. Using h or k (see below)
"random" random blocking. Using k (see below)
"user" user defined groups. Using groups (see below)
Further seven parameters can be specified. "k" determnined the number of groups, "h" the strength + of the correlation used to group parameter and "groups" is used for user defined groups. + "groups" is a vector containing the group number for each parameter. E.g. for three parameters + with the first two in one group, "groups" would be c(1,1,2). + Further pSel and pGroup can be used to influence the choice of groups. In the sampling process + a number of groups is randomly drawn and updated. pSel is a vector containing relative probabilities + for an update of the respective number of groups. E.g. for always updating only one group pSel = 1. + For updating one or two groups with the same probability pSel = c(1,1). By default all numbers + have the same probability. + The same principle is used in pGroup. Here the user can influence the probability of each group + to be updated. By default all groups have the same probability. + Finally "groupStart" defines the starting point of the groupUpdate and "groupIntervall" the intervall + in which the groups are evaluated.
+ +Braak, Cajo JF Ter. "A Markov Chain Monte Carlo version of the genetic algorithm Differential Evolution: easy Bayesian computing for real parameter spaces." Statistics and Computing 16.3 (2006): 239-249.
+ +Differential-Evolution MCMC zs
+ + +DEzs(bayesianSetup, settings = list(iterations = 10000, Z = NULL, startValue = + NULL, pSnooker = 0.1, burnin = 0, thin = 1, f = 2.38, eps = 0, parallel = + NULL, pGamma1 = 0.1, eps.mult = 0.2, eps.add = 0, consoleUpdates = 100, + zUpdateFrequency = 1, currentChain = 1, blockUpdate = list("none", k = NULL, h + = NULL, pSel = NULL, pGroup = NULL, groupStart = 1000, groupIntervall = 1000), + message = TRUE))+ +
bayesianSetup | +a BayesianSetup with the posterior density function to be sampled from |
+
---|---|
settings | +list with parameter settings |
+
startValue | +(optional) eiter a matrix with start population, a number to define the number of chains that are run or a function that samples a starting population. |
+
Z | +starting Z population |
+
iterations | +iterations to run |
+
pSnooker | +probability of Snooker update |
+
burnin | +number of iterations treated as burn-in. These iterations are not recorded in the chain. |
+
thin | +thinning parameter. Determines the interval in which values are recorded. |
+
eps | +small number to avoid singularity |
+
f | +scaling factor gamma |
+
parallel | +logical, determines weather parallel computing should be attempted (see details) |
+
pGamma1 | +probability determining the frequency with which the scaling is set to 1 (allows jumps between modes) |
+
eps.mult | +random term (multiplicative error) |
+
eps.add | +random term |
+
blockUpdate | +list determining whether parameters should be updated in blocks. For possible settings see Details. |
+
message | +logical determines whether the sampler's progress should be printed |
+
For parallel computing FUN needs to take a matrix of proposals
+For blockUpdate the first element in the list determines the type of blocking. +Possible choices are
"none" (default), no blocking of parameters
"correlation" blocking based on correlation of parameters. Using h or k (see below)
"random" random blocking. Using k (see below)
"user" user defined groups. Using groups (see below)
Further seven parameters can be specified. "k" determnined the number of groups, "h" the strength + of the correlation used to group parameter and "groups" is used for user defined groups. + "groups" is a vector containing the group number for each parameter. E.g. for three parameters + with the first two in one group, "groups" would be c(1,1,2). + Further pSel and pGroup can be used to influence the choice of groups. In the sampling process + a number of groups is randomly drawn and updated. pSel is a vector containing relative probabilities + for an update of the respective number of groups. E.g. for always updating only one group pSel = 1. + For updating one or two groups with the same probability pSel = c(1,1). By default all numbers + have the same probability. + The same principle is used in pGroup. Here the user can influence the probability of each group + to be updated. By default all groups have the same probability. + Finally "groupStart" defines the starting point of the groupUpdate and "groupIntervall" the intervall + in which the groups are evaluated.
+ +ter Braak C. J. F., and Vrugt J. A. (2008). Differential Evolution Markov Chain with snooker updater and fewer chains. Statistics and Computing http://dx.doi.org/10.1007/s11222-008-9104-9
+ +Deviance information criterion
+ + +DIC(sampler, ...)+ +
sampler | +An object of class bayesianOutput (mcmcSampler, smcSampler, or mcmcList) |
+
---|---|
... | +further arguments passed to |
+
Output: + list with the following elements: + DIC : Deviance Information Criterion + IC : Bayesian Predictive Information Criterion + pD : Effective number of parameters (pD = Dbar - Dhat) + pV : Effective number of parameters (pV = var(D)/2) + Dbar : Expected value of the deviance over the posterior + Dhat : Deviance at the mean posterior estimate
+ +Spiegelhalter, D. J.; Best, N. G.; Carlin, B. P. & van der Linde, A. (2002) Bayesian measures of model complexity and fit. J. Roy. Stat. Soc. B, 64, 583-639. +Gelman, A.; Hwang, J. & Vehtari, A. (2014) Understanding predictive information criteria for Bayesian models. Statistics and Computing, Springer US, 24, 997-1016-.
+ +The Delayed Rejection Algorithm (Tierney and Mira, 1999)
+ + +DR(startValue = NULL, iterations = 10000, nBI = 0, parmin = NULL, + parmax = NULL, f1 = 1, f2 = 0.5, FUN)+ +
startValue | +vector with the start values for the algorithm. Can be NULL if FUN is of class BayesianSetup. In this case startValues are sampled from the prior. |
+
---|---|
iterations | +iterations to run |
+
nBI | +number of burnin |
+
parmin | +minimum values for the parameter vector or NULL if FUN is of class BayesianSetup |
+
parmax | +maximum values for the parameter vector or NULL if FUN is of class BayesianSetup |
+
f1 | +scaling factor for first proposal |
+
f2 | +scaling factor for second proposal |
+
FUN | +function to be sampled from or object of class bayesianSetup |
+
Tierney, Luke, and Antonietta Mira. "Some adaptive Monte Carlo methods for Bayesian inference." Statistics in medicine 18.1718 (1999): 2507-2515.
+ + +The Delayed Rejection Adaptive Metropolis Algorithm (Haario et al. 2001)
+ + +DRAM(startValue = NULL, iterations = 10000, nBI = 0, parmin = NULL, + parmax = NULL, FUN, f = 1, eps = 0)+ +
startValue | +vector with the start values for the algorithm. Can be NULL if FUN is of class BayesianSetup. In this case startValues are sampled from the prior. |
+
---|---|
iterations | +iterations to run |
+
nBI | +number of burnin |
+
parmin | +minimum values for the parameter vector or NULL if FUN is of class BayesianSetup |
+
parmax | +maximum values for the parameter vector or NULL if FUN is of class BayesianSetup |
+
FUN | +function to be sampled from |
+
f | +scaling factor |
+
eps | +small number to avoid singularity or object of class bayesianSetup |
+
Haario, Heikki, Eero Saksman, and Johanna Tamminen. "An adaptive Metropolis algorithm." Bernoulli (2001): 223-242.
+ + +DREAM
+ + +DREAM(bayesianSetup, settings = list(iterations = 10000, nCR = 3, gamma = + NULL, eps = 0, e = 0.05, pCRupdate = TRUE, updateInterval = 10, burnin = 0, + thin = 1, adaptation = 0.2, DEpairs = 2, consoleUpdates = 10, startValue = + NULL, currentChain = 1, message = TRUE))+ +
bayesianSetup | +Object of class 'bayesianSetup' or 'bayesianOuput'. |
+
---|---|
settings | +list with parameter values |
+
iterations | +Number of model evaluations |
+
nCR | +parameter determining the number of cross-over proposals. If nCR = 1 all parameters are updated jointly. |
+
updateInterval | +determining the intervall for the pCR update |
+
gamma | +Kurtosis parameter Bayesian Inference Scheme. Default is 0. |
+
eps | +Ergodicity term. Default to 0. |
+
e | +Ergodicity term. Default to 5e-2. |
+
pCRupdate | +Update of crossover probabilities, default is TRUE |
+
burnin | +number of iterations treated as burn-in. These iterations are not recorded in the chain. |
+
thin | +thin thinning parameter. Determines the interval in which values are recorded. |
+
adaptation | +Number or percentage of samples that are used for the adaptation in DREAM (see Details). |
+
DEpairs | +Number of pairs used to generate proposal |
+
startValue | +eiter a matrix containing the start values (see details), an integer to define the number of chains that are run, a function to sample the start values or NUll, in which case the values are sampled from the prior. |
+
consoleUpdates | +Intervall in which the sampling progress is printed to the console |
+
message | +logical determines whether the sampler's progress should be printed |
+
mcmc.object containing the following elements: chains, X, pCR
+ +Insted of a bayesianSetup, the function can take the output of a previous run to restart the sampler +from the last iteration. Due to the sampler's internal structure you can only use the output +of DREAM. +If you provide a matrix with start values the number of rows determines the number of chains that are run. +The number of coloumns must be equivalent to the number of parameters in your bayesianSetup. +There are several small differences in the algorithm presented here compared to the original paper by Vrugt et al. (2009). Mainly +the algorithm implemented here does not have an automatic stopping criterion. Hence, it will +always run the number of iterations specified by the user. Also, convergence is not +monitored and left to the user. This can easily be done with coda::gelman.diag(chain). +Further the proposed delayed rejectio step in Vrugt et al. (2009) is not implemented here.
+During the adaptation phase DREAM is running two mechanisms to enhance the sampler's efficiency. +First the disribution of crossover values is tuned to favor large jumps in the parameter space. +The crossover probabilities determine how many parameters are updated simultaneously. +Second outlier chains are replanced as they can largely deteriorate the sampler's performance. +However, these steps destroy the detailed balance of the chain. Consequently these parts of the chain +should be discarded when summarizing posterior moments. This can be done automatically during the +sampling process (i.e. burnin > adaptation) or subsequently by the user. We chose to distinguish between +the burnin and adaptation phase to allow the user more flexibility in the sampler's settings.
+ +Vrugt, Jasper A., et al. "Accelerating Markov chain Monte Carlo simulation by differential evolution with self-adaptive randomized subspace sampling." International Journal of Nonlinear Sciences and Numerical Simulation 10.3 (2009): 273-290.
+ +DREAMzs
+ + +DREAMzs(bayesianSetup, settings = list(iterations = 10000, nCR = 3, gamma = + NULL, eps = 0, e = 0.05, pCRupdate = FALSE, updateInterval = 10, burnin = 0, + thin = 1, adaptation = 0.2, parallel = NULL, Z = NULL, ZupdateFrequency = 10, + pSnooker = 0.1, DEpairs = 2, consoleUpdates = 10, startValue = NULL, + currentChain = 1, message = FALSE))+ +
bayesianSetup | +Object of class 'bayesianSetup' or 'bayesianOuput'. |
+
---|---|
settings | +list with parameter values |
+
iterations | +Number of model evaluations |
+
nCR | +parameter determining the number of cross-over proposals. If nCR = 1 all parameters are updated jointly. |
+
updateInterval | +determining the intervall for the pCR (crossover probabilities) update |
+
gamma | +Kurtosis parameter Bayesian Inference Scheme. Default is 0. |
+
eps | +Ergodicity term. Default to 0. |
+
e | +Ergodicity term. Default to 5e-2. |
+
pCRupdate | +Update of crossover probabilities, default is TRUE |
+
burnin | +number of iterations treated as burn-in. These iterations are not recorded in the chain. |
+
thin | +thin thinning parameter. Determines the interval in which values are recorded. |
+
adaptation | +Number or percentage of samples that are used for the adaptation in DREAM (see Details) |
+
DEpairs | +Number of pairs used to generate proposal |
+
ZupdateFrequency | +frequency to update Z matrix |
+
pSnooker | +probability of snooker update |
+
Z | +starting matrix for Z |
+
startValue | +eiter a matrix containing the start values (see details), an integer to define the number of chains that are run, a function to sample the start values or NUll, in which case the values are sampled from the prior. |
+
consoleUpdates | +Intervall in which the sampling progress is printed to the console |
+
message | +logical determines whether the sampler's progress should be printed |
+
mcmc.object containing the following elements: chains, X, pCR, Z
+ +Insted of a bayesianSetup, the function can take the output of a previous run to restart the sampler +from the last iteration. Due to the sampler's internal structure you can only use the output +of DREAMzs. +If you provide a matrix with start values the number of rows detemines the number of chains that are run. +The number of coloumns must be equivalent to the number of parameters in your bayesianSetup. +There are several small differences in the algorithm presented here compared to the original paper by Vrugt et al. (2009). Mainly +the algorithm implemented here does not have an automatic stopping criterion. Hence, it will +always run the number of iterations specified by the user. Also, convergence is not +monitored and left to the user. This can easily be done with coda::gelman.diag(chain). +Further the proposed delayed rejectio step in Vrugt et al. (2009) is not implemented here. +During the adaptation phase DREAM is running two mechanisms to enhance the sampler's efficiency. +First the disribution of crossover values is tuned to favor large jumps in the parameter space. +The crossover probabilities determine how many parameters are updated simultaneously. +Second outlier chains are replanced as they can largely deteriorate the sampler's performance. +However, these steps destroy the detailed balance of the chain. Consequently these parts of the chain +should be discarded when summarizing posterior moments. This can be done automatically during the +sampling process (i.e. burnin > adaptation) or subsequently by the user. We chose to distinguish between +the burnin and adaptation phase to allow the user more flexibility in the sampler's settings.
+ +Vrugt, Jasper A., et al. "Accelerating Markov chain Monte Carlo simulation by differential evolution with self-adaptive randomized subspace sampling." International Journal of Nonlinear Sciences and Numerical Simulation 10.3 (2009): 273-290.
+ter Braak C. J. F., and Vrugt J. A. (2008). Differential Evolution Markov Chain with snooker updater and fewer chains. Statistics and Computing http://dx.doi.org/10.1007/s11222-008-9104-9
+ +Standard GOF metrics +Startvalues for sampling with nrChains > 1 : if you want to provide different start values for the different chains, provide a list
+ + +GOF(observed, predicted, plot = F)+ +
observed | +observed values |
+
---|---|
predicted | +predicted values |
+
plot | +should a plot be created |
+
A list with the following entries: rmse = root mean squared error, mae = mean absolute error, slope = slope of regression, offset = intercept of regression, R2 = R2 of regression
+ +The function considers predicted ~ observed and calculates
+1) rmse = root mean squared error +2) mae = mean absolute errorr +3) a linear regression with slope, intercept and coefficient of determination R2
+For the linear regression, the x axis is centered, meaning that the intercept is the difference between observed / predicted for the MEAN predicted value. This setting avoids a correlation between slope and intercept (that the intercept is != 0 as soon as the slope is !=0)
+ + +++x = runif(500,-1,1) +y = 0.2 + 0.9 *x + rnorm(500, sd = 0.5) + +summary(lm(y ~ x))#> +#> Call: +#> lm(formula = y ~ x) +#> +#> Residuals: +#> Min 1Q Median 3Q Max +#> -1.44986 -0.30138 -0.02512 0.33049 1.66535 +#> +#> Coefficients: +#> Estimate Std. Error t value Pr(>|t|) +#> (Intercept) 0.20840 0.02258 9.23 <2e-16 *** +#> x 0.88446 0.03645 24.27 <2e-16 *** +#> --- +#> Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 +#> +#> Residual standard error: 0.5041 on 498 degrees of freedom +#> Multiple R-squared: 0.5418, Adjusted R-squared: 0.5409 +#> F-statistic: 588.9 on 1 and 498 DF, p-value: < 2.2e-16 +#>+GOF(x,y)#> $rmse +#> [1] 0.5506948 +#> +#> $mae +#> [1] 0.4308861 +#> +#> $slope +#> [1] 0.8844614 +#> +#> $offset +#> [1] 0.17963 +#> +#> $R2 +#> [1] 0.5418068 +#>+GOF(x,y, plot = TRUE)#> $rmse +#> [1] 0.5506948 +#> +#> $mae +#> [1] 0.4308861 +#> +#> $slope +#> [1] 0.8844614 +#> +#> $offset +#> [1] 0.17963 +#> +#> $R2 +#> [1] 0.5418068 +#>
The Metropolis Algorithm (Metropolis et al. 1953)
+ + +M(startValue = NULL, iterations = 10000, nBI = 0, parmin = NULL, + parmax = NULL, f = 1, FUN, consoleUpdates = 1000)+ +
startValue | +vector with the start values for the algorithm. Can be NULL if FUN is of class BayesianSetup. In this case startValues are sampled from the prior. |
+
---|---|
iterations | +iterations to run |
+
nBI | +number of burnin |
+
parmin | +minimum values for the parameter vector or NULL if FUN is of class BayesianSetup |
+
parmax | +maximum values for the parameter vector or NULL if FUN is of class BayesianSetup |
+
f | +scaling factor |
+
FUN | +function to be sampled from or object of class bayesianSetup |
+
consoleUpdates | +interger, determines the frequency with which sampler progress is printed to the console |
+
Metropolis, Nicholas, et al. "Equation of state calculations by fast computing machines." The journal of chemical physics 21.6 (1953): 1087-1092.
+ + +calculates the Maxiumum APosteriori value (MAP)
+ + +MAP(bayesianOutput, ...)+ +
bayesianOutput | +an object of class BayesianOutput (mcmcSampler, smcSampler, or mcmcList) |
+
---|---|
... | +optional values to be passed on the the getSample function |
+
Currently, this function simply returns the parameter combination with the highest posterior in the chain. A more refined option would be to take the MCMC sample and do additional calculations, e.g. use an optimizer, a kerne delnsity estimator, or some other tool to search / interpolate around the best value in the chain
+ +Creates a Metropolis-type MCMC with options for covariance adaptatin, delayed rejection, Metropolis-within-Gibbs, and tempering
+ + +Metropolis(bayesianSetup, settings = list(startValue = NULL, optimize = T, + proposalGenerator = NULL, consoleUpdates = 100, burnin = 0, thin = 1, parallel + = NULL, adapt = T, adaptationInterval = 500, adaptationNotBefore = 3000, + DRlevels = 1, proposalScaling = NULL, adaptationDepth = NULL, + temperingFunction = NULL, gibbsProbabilities = NULL, message = TRUE))+ +
bayesianSetup | +either an object of class bayesianSetup created by |
+
---|---|
settings | +a list of settings - possible options follow below |
+
startValue | +startValue for the MCMC and optimization (if optimize = T). If not provided, the sampler will attempt to obtain the startValue from the bayesianSetup |
+
optimize | +logical, determines whether an optimization for start values and proposal function should be run before starting the sampling |
+
proposalGenerator | +optional proposalgenerator object (see |
+
proposalScaling | +additional scaling parameter for the proposals that controls the different scales of the proposals after delayed rejection (typical, after a rejection, one would want to try a smaller scale). Needs to be as long as DRlevels. Defaults to 0.5^(- 0:(mcmcSampler$settings$DRlevels -1) |
+
burnin | +number of iterations treated as burn-in. These iterations are not recorded in the chain. |
+
thin | +thinning parameter. Determines the interval in which values are recorded. |
+
consoleUpdates | +integer, determines the frequency with which sampler progress is printed to the console |
+
adapt | +logical, determines wheter an adaptive algorithm should be implemented. Default is TRUE. |
+
adaptationInterval | +integer, determines the interval of the adaption if adapt = TRUE. |
+
adaptationNotBefore | +integer, determines the start value for the adaption if adapt = TRUE. |
+
DRlevels | +integer, determines the number of levels for a delayed rejection sampler. Default is 1, which means no delayed rejection is used. |
+
temperingFunction | +function to implement simulated tempering in the algorithm. The function describes how the acceptance rate will be influenced in the course of the iterations. |
+
gibbsProbabilities | +vector that defines the relative probabilities of the number of parameters to be changes simultaniously. |
+
message | +logical determines whether the sampler's progress should be printed |
+
The 'Metropolis' function is the main function for all Metropolis based samplers in this package. To call the derivatives from the basic Metropolis-Hastings MCMC, you can either use the corresponding function (e.g. AM
for an adaptive Metropolis sampler) or use the parameters to adapt the basic Metropolis-Hastings. The advantage of the latter case is that you can easily combine different properties (e.g. adapive sampling and delayed rejection sampling) without changing the function.
Haario, H., E. Saksman, and J. Tamminen (2001). An adaptive metropolis algorithm. Bernoulli , 223-242.
+Haario, Heikki, et al. "DRAM: efficient adaptive MCMC." Statistics and Computing 16.4 (2006): 339-354.
+Hastings, W. K. (1970). Monte carlo sampling methods using markov chains and their applications. Biometrika 57 (1), 97-109.
+Green, Peter J., and Antonietta Mira. "Delayed rejection in reversible jump Metropolis-Hastings." Biometrika (2001): 1035-1053.
+Metropolis, N., A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller (1953). Equation of state calculations by fast computing machines. The journal of chemical physics 21 (6), 1087 - 1092.
+ + +++# Running the metropolis via the runMCMC with a proposal covariance generated from the prior +# (can be useful for complicated priors) + +ll = function(x) sum(dnorm(x, log = TRUE)) +setup = createBayesianSetup(ll, lower = c(-10,-10), upper = c(10,10)) + +samples = setup$prior$sampler(1000) + +generator = createProposalGenerator(diag(1, setup$numPars)) +generator = updateProposalGenerator(generator, samples, manualScaleAdjustment = 1, message = TRUE)#> Proposalgenerator settings changedcovariance set to: +#> [,1] [,2] +#> [1,] 34.8160170 -0.5583527 +#> [2,] -0.5583527 32.8871576 +#> covarianceDecomp set to: +#> [,1] [,2] +#> [1,] 5.90051 -0.09462788 +#> [2,] 0.00000 5.73395180 +#> gibbsProbabilities set to: +#> NULL +#> gibbsWeights set to: +#> NULL +#> otherDistribution set to: +#> NULL +#> otherDistributionLocation set to: +#> NULL+settings = list(proposalGenerator = generator, optimize = FALSE, iterations = 500) + +out = runMCMC(bayesianSetup = setup, sampler = "Metropolis", settings = settings)#> Running Metropolis-MCMC, chain 1 iteration 100 of 500 . Current logp: -12.60143 Please wait! Running Metropolis-MCMC, chain 1 iteration 200 of 500 . Current logp: -8.423186 Please wait! Running Metropolis-MCMC, chain 1 iteration 300 of 500 . Current logp: -8.174177 Please wait! Running Metropolis-MCMC, chain 1 iteration 400 of 500 . Current logp: -8.532211 Please wait! Running Metropolis-MCMC, chain 1 iteration 500 of 500 . Current logp: -8.532211 Please wait!#>
T-walk MCMC
+ + +Twalk(bayesianSetup, settings = list(iterations = 10000, at = 6, aw = 1.5, pn1 + = NULL, Ptrav = 0.4918, Pwalk = 0.4918, Pblow = 0.0082, burnin = 0, thin = 1, + startValue = NULL, consoleUpdates = 100, message = TRUE))+ +
bayesianSetup | +Object of class 'bayesianSetup' or 'bayesianOuput'. |
+
---|---|
settings | +list with parameter values. |
+
iterations | +Number of model evaluations |
+
at | +"traverse" move proposal parameter. Default to 6 |
+
aw | +"walk" move proposal parameter. Default to 1.5 |
+
pn1 | +Probability determining the number of parameters that are changed |
+
Ptrav | +Move probability of "traverse" moves, default to 0.4918 |
+
Pwalk | +Move probability of "walk" moves, default to 0.4918 |
+
Pblow | +Move probability of "traverse" moves, default to 0.0082 |
+
burnin | +number of iterations treated as burn-in. These iterations are not recorded in the chain. |
+
thin | +thinning parameter. Determines the interval in which values are recorded. |
+
startValue | +Matrix with start values |
+
consoleUpdates | +Intervall in which the sampling progress is printed to the console |
+
message | +logical determines whether the sampler's progress should be printed |
+
Object of class bayesianOutput.
+ +The probability of "hop" moves is 1 minus the sum of all other probabilities.
+ +Christen, J. Andres, and Colin Fox. "A general purpose sampling algorithm for continuous distributions (the t-walk)." Bayesian Analysis 5.2 (2010): 263-281.
+ + +A very simple ecosystem model, based on three carbon pools and a basic LUE model
+ + +VSEM(pars = c(KEXT = 0.5, LAR = 1.5, LUE = 0.002, GAMMA = 0.4, tauV = 1440, + tauS = 27370, tauR = 1440, Av = 0.5, Cv = 3, Cs = 15, Cr = 3), PAR, + C = TRUE)+ +
pars | +a parameter vector with parameters and initial states |
+
---|---|
PAR | +Forcing, photosynthetically active radiation (PAR) MJ /m2 /day |
+
C | +switch to choose whether to use the C or R version of the model. C is much faster. |
+
a matrix with colums NEE, CV, CR and CS units and explanations see details
+ +This Very Simple Ecosystem Model (VSEM) is a 'toy' model designed to be very simple but yet bear some resemblance to deterministic processed based ecosystem models (PBMs) that are commonly used in forest modelling.
+The model determines the accumulation of carbon in the plant and soil from the growth of the plant via photosynthesis and senescence to the soil which respires carbon back to the atmosphere.
+The model calculates Gross Primary Productivity (GPP) using a very simple light-use efficiency (LUE) formulation multiplied by light interception. Light interception is calculated via Beer's law with a constant light extinction coefficient operating on Leaf Area Index (LAI).
+A parameter (GAMMA) determines the fraction of GPP that is autotrophic respiration. The Net Primary Productivity (NPP) is then allocated to above and below-ground vegetation via a fixed allocation fraction. Carbon is lost from the plant pools to a single soil pool via fixed turnover rates. Heterotropic respiration in the soil is determined via a soil turnover rate.
+The model equations are
+-- Photosynthesis + $$LAI = LAR*Cv$$ +$$GPP = PAR * LUE * (1 - \exp^{(-KEXT * LAI)})$$ +$$NPP = (1-GAMMA) * GPP$$
+-- State equations +$$dCv/dt = Av * NPP - Cv/tauV$$ +$$dCr/dt = (1.0-Av) * NPP - Cr/tauR$$ +$$dCs/dt = Cr/tauR + Cv/tauV - Cs/tauS$$
+The model time-step is daily.
+-- VSEM inputs:
+PAR Photosynthetically active radiation (PAR) MJ /m2 /day
+-- VSEM parameters:
+KEXT Light extinction coefficient m2 ground area / m2 leaf area
+LAR Leaf area ratio m2 leaf area / kg aboveground vegetation
+LUE Light-Use Efficiency (kg C MJ-1 PAR)
+GAMMA Autotrophic respiration as a fraction of GPP
+tauV Longevity of aboveground vegetation days
+tauR Longevity of belowground vegetation days
+tauS Residence time of soil organic matter d
+-- VSEM states:
+Cv Above-ground vegetation pool kg C / m2
+Cr Below-ground vegetation pool kg C / m2
+Cs Carbon in organic matter kg C / m2
+-- VSEM fluxes:
+G Gross Primary Productivity kg C /m2 /day
+NPP Net Primary Productivity kg C /m2 /day
+NEE Net Ecosystem Exchange kg C /m2 /day
+ +VSEMgetDefaults
, VSEMcreatePAR
, , VSEMcreateLikelihood
++ +## This example shows how to run and calibrate the VSEM model + +library(BayesianTools) + +# Create input data for the model +PAR <- VSEMcreatePAR(1:1000) +plot(PAR, main = "PAR (driving the model)", xlab = "Day")+# load reference parameter definition (upper, lower prior) +refPars <- VSEMgetDefaults() +# this adds one additional parameter for the likelihood standard deviation (see below) +refPars[12,] <- c(2, 0.1, 4) +rownames(refPars)[12] <- "error-sd" +head(refPars)#> best lower upper +#> KEXT 0.500 2e-01 1e+00 +#> LAR 1.500 2e-01 3e+00 +#> LUE 0.002 5e-04 4e-03 +#> GAMMA 0.400 2e-01 6e-01 +#> tauV 1440.000 5e+02 3e+03 +#> tauS 27370.000 4e+03 5e+04+# create some simulated test data +# generally recommended to start with simulated data before moving to real data +referenceData <- VSEM(refPars$best[1:11], PAR) # model predictions with reference parameters +referenceData[,1] = 1000 * referenceData[,1] +# this adds the error - needs to conform to the error definition in the likelihood +obs <- referenceData + rnorm(length(referenceData), sd = refPars$best[12]) +oldpar <- par(mfrow = c(2,2)) +for (i in 1:4) plotTimeSeries(observed = obs[,i], + predicted = referenceData[,i], main = colnames(referenceData)[i])+# Best to program in a way that we can choose easily which parameters to calibrate +parSel = c(1:6, 12) + +# here is the likelihood +likelihood <- function(par, sum = TRUE){ + # set parameters that are not calibrated on default values + x = refPars$best + x[parSel] = par + predicted <- VSEM(x[1:11], PAR) # replace here VSEM with your model + predicted[,1] = 1000 * predicted[,1] # this is just rescaling + diff <- c(predicted[,1:4] - obs[,1:4]) # difference betweeno observed and predicted + # univariate normal likelihood. Note that there is a parameter involved here that is fit + llValues <- dnorm(diff, sd = x[12], log = TRUE) + if (sum == FALSE) return(llValues) + else return(sum(llValues)) +} + +# optional, you can also directly provide lower, upper in the createBayesianSetup, see help +prior <- createUniformPrior(lower = refPars$lower[parSel], + upper = refPars$upper[parSel], best = refPars$best[parSel]) + +bayesianSetup <- createBayesianSetup(likelihood, prior, names = rownames(refPars)[parSel]) + +# settings for the sampler, iterations should be increased for real applicatoin +settings <- list(iterations = 2000, nrChains = 2) + +out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "DEzs", settings = settings)#> Running DEzs-MCMC, chain 1 iteration 300 of 2001 . Current logp -8624.918 -8592.876 -8588.121 . Please wait! Running DEzs-MCMC, chain 1 iteration 600 of 2001 . Current logp -8532.735 -8537.345 -8515.801 . Please wait! Running DEzs-MCMC, chain 1 iteration 900 of 2001 . Current logp -8519.858 -8525.309 -8515.801 . Please wait! Running DEzs-MCMC, chain 1 iteration 1200 of 2001 . Current logp -8506.742 -8506.346 -8515.801 . Please wait! Running DEzs-MCMC, chain 1 iteration 1500 of 2001 . Current logp -8504.459 -8506.262 -8501.071 . Please wait! Running DEzs-MCMC, chain 1 iteration 1800 of 2001 . Current logp -8505.372 -8499.451 -8500.069 . Please wait! Running DEzs-MCMC, chain 1 iteration 2001 of 2001 . Current logp -8503.789 -8499.372 -8498.692 . Please wait!#>#> Running DEzs-MCMC, chain 2 iteration 300 of 2001 . Current logp -8507.99 -8534.268 -8518.524 . Please wait! Running DEzs-MCMC, chain 2 iteration 600 of 2001 . Current logp -8504.173 -8511.363 -8503.835 . Please wait! Running DEzs-MCMC, chain 2 iteration 900 of 2001 . Current logp -8497.267 -8502.677 -8497.848 . Please wait! Running DEzs-MCMC, chain 2 iteration 1200 of 2001 . Current logp -8497.805 -8502.279 -8497.848 . Please wait! Running DEzs-MCMC, chain 2 iteration 1500 of 2001 . Current logp -8496.634 -8497.897 -8497.356 . Please wait! Running DEzs-MCMC, chain 2 iteration 1800 of 2001 . Current logp -8499.464 -8498.916 -8497.375 . Please wait! Running DEzs-MCMC, chain 2 iteration 2001 of 2001 . Current logp -8499.549 -8500.452 -8499.336 . Please wait!#>+# NOT RUN { +plot(out) +summary(out) +marginalPlot(out, scale = T) +gelmanDiagnostics(out) # should be below 1.05 for all parameters to demonstrate convergence + +# Posterior predictive simulations + +# Create a prediction function +createPredictions <- function(par){ + # set the parameters that are not calibrated on default values + x = refPars$best + x[parSel] = par + predicted <- VSEM(x[1:11], PAR) # replace here VSEM with your model + return(predicted[,1] * 1000) +} + +# Create an error function +createError <- function(mean, par){ + return(rnorm(length(mean), mean = mean, sd = par[7])) +} + +# plot prior predictive distribution and prior predictive simulations +plotTimeSeriesResults(sampler = out, model = createPredictions, observed = obs[,1], + error = createError, prior = TRUE, main = "Prior predictive") + +# plot posterior predictive distribution and posterior predictive simulations +plotTimeSeriesResults(sampler = out, model = createPredictions, observed = obs[,1], + error = createError, main = "Posterior predictive") + + +######################################################## +# Demonstrating the updating of the prior from old posterior +# Note that it is usually more exact to rerun the MCMC +# with all (old and new) data, instead of updating the prior +# because likely some information is lost when approximating the +# Posterior by a multivariate normal + +settings <- list(iterations = 5000, nrChains = 2) + +out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "DEzs", settings = settings) + +plot(out) +correlationPlot(out, start = 1000) + +newPrior = createPriorDensity(out, method = "multivariate", + eps = 1e-10, + lower = refPars$lower[parSel], + upper = refPars$upper[parSel], start= 1000) + +bayesianSetup <- createBayesianSetup(likelihood = likelihood, + prior = newPrior, + names = rownames(refPars)[parSel] ) + +# check boundaries are correct set +bayesianSetup$prior$sampler() < refPars$lower[parSel] +bayesianSetup$prior$sampler() > refPars$upper[parSel] + +# check prior looks similar to posterior +x = bayesianSetup$prior$sampler(2000) +correlationPlot(x, thin = F) + +out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "DEzs", settings = settings) + +plot(out) +correlationPlot(out) + +plotTimeSeriesResults(sampler = out, + model = createPredictions, + observed = obs[,1], + error = createError, + prior = F, main = "Posterior predictive") + +plotTimeSeriesResults(sampler = out, + model = createPredictions, + observed = obs[,1], + error = createError, + prior = T, main = "Prior predictive") + + + + +# }+par(oldpar)
Create an example dataset, and from that a likelihood or posterior for the VSEM model
+ + +VSEMcreateLikelihood(likelihoodOnly = F, plot = F, selection = c(1:6, 12))+ +
likelihoodOnly | +switch to devide whether to create only a likelihood, or a full bayesianSetup with uniform priors. |
+
---|---|
plot | +switch to decide whether data should be plotted |
+
selection | +vector containing the indices of the selected parameters |
+
The purpose of this function is to be able to conveniently create a likelihood for the VSEM model for demonstration purposes. The function creates example data --> likelihood --> BayesianSetup, where the latter is the
+ + +calculates the WAIC
+ + +WAIC(bayesianOutput, numSamples = 1000, ...)+ +
bayesianOutput | +an object of class BayesianOutput |
+
---|---|
numSamples | +the number of samples to calculate the WAIC |
+
... | +optional values to be passed on the the getSample function |
+
The WAIC is constructed as
+The lppd (log pointwise predictive density), defined in Gelman et al., 2013, eq. 4 as
+The value of can be calculated in two ways, the method used is determined by the
+method
argument.
Method 1 is defined as, + Method 2 is defined as, + where is the sample variance.
+ +The function requires that the likelihood passed on to BayesianSetup contains the option sum = T/F, with defaul F. If set to true, the likelihood for each data point must be returned.
+ +Gelman, Andrew and Jessica Hwang and Aki Vehtari (2013), "Understanding Predictive Information Criteria for Bayesian Models," http://www.stat.columbia.edu/~gelman/research/unpublished/waic_understand_final.pdf.
+Watanabe, S. (2010). "Asymptotic Equivalence of Bayes Cross Validation and Widely Applicable Information Criterion in Singular Learning Theory", Journal of Machine Learning Research, http://www.jmlr.org/papers/v11/watanabe10a.html.
+ ++bayesianSetup <- createBayesianSetup(likelihood = testDensityNormal, + prior = createUniformPrior(lower = rep(-10,2), + upper = rep(10,2))) + +# likelihood density needs to have option sum = FALSE + +testDensityNormal(c(1,1,1), sum = FALSE)#> [1] -1.418939 -1.418939 -1.418939bayesianSetup$likelihood$density(c(1,1,1), sum = FALSE)#> [1] -1.418939 -1.418939 -1.418939bayesianSetup$likelihood$density(matrix(rep(1,9), ncol = 3), sum = FALSE)#> [,1] [,2] [,3] +#> [1,] -1.418939 -1.418939 -1.418939 +#> [2,] -1.418939 -1.418939 -1.418939 +#> [3,] -1.418939 -1.418939 -1.418939#> Running DEzs-MCMC, chain 1 iteration 300 of 10002 . Current logp -8.435686 -7.89902 -8.224242 . Please wait! Running DEzs-MCMC, chain 1 iteration 600 of 10002 . Current logp -9.623622 -8.0001 -9.358736 . Please wait! Running DEzs-MCMC, chain 1 iteration 900 of 10002 . Current logp -8.210087 -9.184124 -7.869296 . Please wait! Running DEzs-MCMC, chain 1 iteration 1200 of 10002 . Current logp -8.680172 -8.612244 -8.269293 . Please wait! Running DEzs-MCMC, chain 1 iteration 1500 of 10002 . Current logp -9.730421 -8.775348 -8.081375 . Please wait! Running DEzs-MCMC, chain 1 iteration 1800 of 10002 . Current logp -8.270908 -8.957659 -10.43252 . Please wait! Running DEzs-MCMC, chain 1 iteration 2100 of 10002 . Current logp -7.984955 -8.251509 -8.458087 . Please wait! Running DEzs-MCMC, chain 1 iteration 2400 of 10002 . Current logp -7.882363 -7.898614 -9.155126 . Please wait! Running DEzs-MCMC, chain 1 iteration 2700 of 10002 . Current logp -9.426532 -8.327137 -8.185499 . Please wait! Running DEzs-MCMC, chain 1 iteration 3000 of 10002 . Current logp -9.454616 -9.183276 -9.804906 . Please wait! Running DEzs-MCMC, chain 1 iteration 3300 of 10002 . Current logp -9.85258 -8.817787 -8.349674 . Please wait! Running DEzs-MCMC, chain 1 iteration 3600 of 10002 . Current logp -7.963722 -8.886764 -8.318415 . Please wait! Running DEzs-MCMC, chain 1 iteration 3900 of 10002 . Current logp -11.09728 -9.129186 -9.112218 . Please wait! Running DEzs-MCMC, chain 1 iteration 4200 of 10002 . Current logp -7.944314 -12.44574 -9.136989 . Please wait! Running DEzs-MCMC, chain 1 iteration 4500 of 10002 . Current logp -8.305425 -7.936811 -7.953184 . Please wait! Running DEzs-MCMC, chain 1 iteration 4800 of 10002 . Current logp -10.98229 -8.645307 -8.081101 . Please wait! Running DEzs-MCMC, chain 1 iteration 5100 of 10002 . Current logp -8.032454 -8.045189 -9.515575 . Please wait! Running DEzs-MCMC, chain 1 iteration 5400 of 10002 . Current logp -8.120442 -10.54898 -9.443533 . Please wait! Running DEzs-MCMC, chain 1 iteration 5700 of 10002 . Current logp -8.049416 -8.717827 -7.862733 . Please wait! Running DEzs-MCMC, chain 1 iteration 6000 of 10002 . Current logp -8.084893 -8.251851 -9.081135 . Please wait! Running DEzs-MCMC, chain 1 iteration 6300 of 10002 . Current logp -8.073244 -7.974431 -8.079902 . Please wait! Running DEzs-MCMC, chain 1 iteration 6600 of 10002 . Current logp -8.218521 -9.22812 -8.772962 . Please wait! Running DEzs-MCMC, chain 1 iteration 6900 of 10002 . Current logp -8.246165 -8.5875 -8.108148 . Please wait! Running DEzs-MCMC, chain 1 iteration 7200 of 10002 . Current logp -10.74085 -7.880688 -10.23808 . Please wait! Running DEzs-MCMC, chain 1 iteration 7500 of 10002 . Current logp -8.916833 -8.728722 -8.486303 . Please wait! Running DEzs-MCMC, chain 1 iteration 7800 of 10002 . Current logp -8.363117 -8.274794 -8.98649 . Please wait! Running DEzs-MCMC, chain 1 iteration 8100 of 10002 . Current logp -9.306803 -8.630766 -9.00387 . Please wait! Running DEzs-MCMC, chain 1 iteration 8400 of 10002 . Current logp -7.83784 -9.938769 -7.96229 . Please wait! Running DEzs-MCMC, chain 1 iteration 8700 of 10002 . Current logp -10.16324 -8.538351 -8.001384 . Please wait! Running DEzs-MCMC, chain 1 iteration 9000 of 10002 . Current logp -9.678691 -8.667391 -10.28152 . Please wait! Running DEzs-MCMC, chain 1 iteration 9300 of 10002 . Current logp -8.111272 -10.70121 -8.247235 . Please wait! Running DEzs-MCMC, chain 1 iteration 9600 of 10002 . Current logp -8.228964 -9.800455 -12.53619 . Please wait! Running DEzs-MCMC, chain 1 iteration 9900 of 10002 . Current logp -8.535171 -8.331085 -7.84911 . Please wait! Running DEzs-MCMC, chain 1 iteration 10002 of 10002 . Current logp -8.401052 -8.317495 -8.432323 . Please wait!#>+WAIC(out)#> $WAIC1 +#> [1] 6.345918 +#> +#> $WAIC2 +#> [1] 7.044639 +#> +#> $lppd +#> [1] -2.55475 +#> +#> $pWAIC1 +#> [1] 0.6182095 +#> +#> $pWAIC2 +#> [1] 0.9675698 +#>
Provides the default settings for the different samplers in runMCMC
+ + +applySettingsDefault(settings = NULL, sampler = "DEzs", check = FALSE)+ +
settings | +optional list with parameters that will be used instead of the defaults |
+
---|---|
sampler | +one of the samplers in |
+
check | +logical determines whether parameters should be checked for consistency |
+
see runMCMC
Calculates the marginal likelihood of a chain via bridge sampling
+ + +bridgesample(chain, nParams, lower = NULL, upper = NULL, ...)+ +
chain | +a single mcmc chain with samples as rows and parameters and posterior density as columns. |
+
---|---|
nParams | +number of parameters |
+
lower | +optional - lower bounds of the prior |
+
upper | +optional - upper bounds of the prior |
+
... | +arguments passed to bridge_sampler |
+
This function uses "bridge_sampler" from the package "bridgesampling".
+ ++means <- c(0, 1, 2) +sds <- c(1, 0.6, 3) + +# log-likelihood +ll <- function (x) { + return(sum(dnorm(x, mean = means, sd = sds, log = TRUE))) +} + +# lower and upper bounds for prior +lb <- rep(-10, 3) +ub <- rep(10, 3) + +# create setup and run MCMC +setup <- createBayesianSetup(likelihood = ll, + lower = lb, + upper = ub) + +out <- runMCMC(bayesianSetup = setup, + settings = list(iterations = 1000), + sampler = "DEzs")#> Running DEzs-MCMC, chain 1 iteration 300 of 1002 . Current logp -14.67601 -13.46839 -16.29187 . Please wait! Running DEzs-MCMC, chain 1 iteration 600 of 1002 . Current logp -14.8849 -14.08858 -15.73934 . Please wait! Running DEzs-MCMC, chain 1 iteration 900 of 1002 . Current logp -15.12236 -14.57723 -12.89714 . Please wait! Running DEzs-MCMC, chain 1 iteration 1002 of 1002 . Current logp -13.07464 -13.45026 -14.50715 . Please wait!#>+# sample from MCMC output with "burn-in" of 25% +sample <- getSample(out$chain, start = 250, numSamples = 500) + +# use bridge sampling to get marginal likelihood +bs_result <- bridgesample(chain = sample, + nParams = out$setup$numPars, + lower = lb, + upper = ub)#> Iteration: 1 +#> Iteration: 2 +#> Iteration: 3 +#> Iteration: 4 +#> Iteration: 5 +#> Iteration: 6 +#> Iteration: 7 +#> Iteration: 8 +#> Iteration: 9 +#> Iteration: 10bs_result#> Bridge sampling estimate of the log marginal likelihood: -9.7342 +#> Estimate obtained in 10 iteration(s) via method "normal".
Function used to assure that an object is of class 'BayesianSetup'. If you pass a function, it is coverted to an object of class 'BayesianSetup' (using createBayesianSetup
) before it is returned.
checkBayesianSetup(bayesianSetup, parallel = F)+ +
bayesianSetup | +either object of class bayesianSetup or a log posterior function |
+
---|---|
parallel | +if bayesianSetup is a function, this will set the parallelization option for the class BayesianSetup that is created internally. If bayesianSetup is already a BayesianSetup, then this will check if parallel = T is requested but not supported by the BayesianSetup. This option is for internal use in the samplers |
+
The recommended option to use this function in the samplers is to have parallel with default NULL in the samplers, so that checkBayesianSetup with a function will create a bayesianSetup without parallelization, while it will do nothing with an existing BayesianSetup. If the user sets parallelization, it will set the approriate parallelization for a function, and check in case of an existing BayesianSetup. The checkBayesianSetup call in the samplers should then be followed by a check for parallel = NULL in sampler, in which case paralell can be set from the BayesianSetup
+ +Function is used to make the plot and diagnostic functions +available for coda::mcmc objects
+ + +convertCoda(sampler, names = NULL, info = NULL, likelihood = NULL)+ +
sampler | +An object of class mcmc or mcmc.list |
+
---|---|
names | +vector giving the parameter names (optional) |
+
info | +matrix (or list with matrices for mcmc.list objects) with three coloumns containing log posterior, log likelihood and log prior of the sampler for each time step (optional; but see Details) |
+
likelihood | +likelihood function used in the sampling (see Details) |
+
The parameter 'likelihood' is optional for most functions but can be needed e.g for
+using the DIC
function.
Also the parameter info is optional for most uses. However for some functions (e.g. MAP
)
+the matrix or single coloumns (e.g. log posterior) are necessary for the diagnostics.
Checks if thin is conistent with nTotalSamples samples and if not corrects it.
+ + +correctThin(nTotalSamples, thin, autoThinFraction = 0.001)+ +
nTotalSamples | +total number of rows/samples |
+
---|---|
thin | +thinning |
+
autoThinFraction | +fraction of the data that will be sampled when thin is set to "auto". E.g. 0.5 means thin will be nTotalSamples * 0.5. The resulting thin value is rounded down to the next integer. |
+
Checks if the thin argument is consistent with the data consisting of nTotalSamples samples/rows and corrects thin if not.
+ + +Flexible function to create correlation density plots
+ + +correlationPlot(mat, density = "smooth", thin = "auto", + method = "pearson", whichParameters = NULL, ...)+ +
mat | +object of class "bayesianOutput" or a matrix or data frame of variables |
+
---|---|
density | +type of plot to do. Either "smooth" (default), "corellipseCor", or "ellipse" |
+
thin | +thinning of the matrix to make things faster. Default is to thin to 5000 |
+
method | +method for calculating correlations. Possible choices are "pearson" (default), "kendall" and "spearman" |
+
whichParameters | +indices of parameters that should be plotted |
+
... | +additional parameters to pass on to the |
+
The code for the correlation density plot originates from Hartig, F.; Dislich, C.; Wiegand, T. & Huth, A. (2014) Technical Note: Approximate Bayesian parameterization of a process-based tropical forest model. Biogeosciences, 11, 1261-1272.
+ +marginalPlot
+ plotTimeSeries
+ tracePlot
+
Creates a standardized collection of prior, likelihood and posterior functions, including error checks etc.
+ + +createBayesianSetup(likelihood, prior = NULL, priorSampler = NULL, + parallel = FALSE, lower = NULL, upper = NULL, best = NULL, + names = NULL, parallelOptions = list(variables = "all", packages = "all", + dlls = NULL), catchDuplicates = FALSE)+ +
likelihood | +log likelihood density function |
+
---|---|
prior | +either a prior class (see |
+
priorSampler | +if a prior density (and not a prior class) is provided to prior, the optional prior sampling function can be provided here |
+
parallel | +parallelization option. Default is F. Other options are T, or "external". See details. |
+
lower | +vector with lower prior limits |
+
upper | +vector with upper prior limits |
+
best | +vector with best values |
+
names | +optional vector with parameter names |
+
parallelOptions | +list containing three lists. First "packages" determines the R packages necessary to run the likelihood function. Second "variables" the objects in the global envirnment needed to run the likelihood function and third "dlls" the DLLs needed to run the likelihood function (see Details and Examples). |
+
catchDuplicates | +Logical, determines whether unique parameter combinations should only be evaluated once. Only used when the likelihood accepts a matrix with parameter as columns. |
+
If lower and upper (and optionally best) but not a prior object are passed, the function will create an uniform prior with the limits lower and upper. If a prior object and lower and upper are passed, the function will ignore lower and upper, and will only use the prior object.
+For parallelization, option T means that an automatic parallelization via R is attempted, or "external", in which case it is assumed that the likelihood is already parallelized. In this case it needs to accept a matrix with parameters as columns. +Further you can specify the packages, objects and DLLs that are exported to the cluster. +By default a copy of your workspace is exported. However, depending on your workspace this can be very inefficient.
+For more details, make sure to read the vignette (run vignette("BayesianTools", package="BayesianTools")
+ +checkBayesianSetup
+ createLikelihood
+ createPrior
+ll <- function(x) sum(dnorm(x, log = TRUE)) + +test <- createBayesianSetup(ll, prior = NULL, priorSampler = NULL, lower = -10, upper = 10) +str(test)#> List of 7 +#> $ prior :List of 6 +#> ..$ density :function (x) +#> .. ..- attr(*, "srcref")=Class 'srcref' atomic [1:8] 59 21 63 3 21 3 59 63 +#> .. .. .. ..- attr(*, "srcfile")=Classes 'srcfilecopy', 'srcfile' <environment: 0x000000002766d2b8> +#> ..$ sampler :function (n = NULL) +#> .. ..- attr(*, "srcref")=Class 'srcref' atomic [1:8] 68 24 76 5 24 5 68 76 +#> .. .. .. ..- attr(*, "srcfile")=Classes 'srcfilecopy', 'srcfile' <environment: 0x000000002766d2b8> +#> ..$ lower : num -10 +#> ..$ upper : num 10 +#> ..$ best : num 0 +#> ..$ originalDensity:function (x) +#> .. ..- attr(*, "srcref")=Class 'srcref' atomic [1:8] 101 14 104 3 14 3 101 104 +#> .. .. .. ..- attr(*, "srcfile")=Classes 'srcfilecopy', 'srcfile' <environment: 0x000000002766d2b8> +#> ..- attr(*, "class")= chr "prior" +#> $ likelihood:List of 3 +#> ..$ density:function (x, ...) +#> .. ..- attr(*, "srcref")=Class 'srcref' atomic [1:8] 45 21 95 3 21 3 45 95 +#> .. .. .. ..- attr(*, "srcfile")=Classes 'srcfilecopy', 'srcfile' <environment: 0x00000000255f6b60> +#> ..$ sampler: NULL +#> ..$ cl : NULL +#> ..- attr(*, "class")= chr "likelihood" +#> $ posterior :List of 1 +#> ..$ density:function (x, returnAll = F) +#> .. ..- attr(*, "srcref")=Class 'srcref' atomic [1:8] 9 16 41 3 16 3 9 41 +#> .. .. .. ..- attr(*, "srcfile")=Classes 'srcfilecopy', 'srcfile' <environment: 0x0000000028eeb968> +#> ..- attr(*, "class")= chr "posterior" +#> $ names : chr "par 1" +#> $ numPars : int 1 +#> $ model : NULL +#> $ parallel : logi FALSE +#> - attr(*, "class")= chr "BayesianSetup"test$prior$density(0)#> [1] -2.995732+test$likelihood$density(c(1,1))#> [1] -2.837877test$likelihood$density(1)#> [1] -1.418939test$posterior$density(1)#> [1] -4.414671test$posterior$density(1, returnAll = TRUE)#> [1] -4.414671 -1.418939 -2.995732+test$likelihood$density(matrix(rep(1,4), nrow = 2))#> [1] -2.837877 -2.837877#test$posterior$density(matrix(rep(1,4), nrow = 2), returnAll = TRUE) +test$likelihood$density(matrix(rep(1,4), nrow = 4))#> [1] -1.418939 -1.418939 -1.418939 -1.418939+ + +## Example of how to use parallelization using the VSEM model +# Note that the parallelization produces overhead and is not always +# speeding things up. In this example, due to the small +# computational cost of the VSEM the parallelization is +# most likely to reduce the speed of the sampler. + +# Creating reference data +PAR <- VSEMcreatePAR(1:1000) +refPars <- VSEMgetDefaults() +refPars[12,] <- c(0.2, 0.001, 1) +rownames(refPars)[12] <- "error-sd" + +referenceData <- VSEM(refPars$best[1:11], PAR) +obs = apply(referenceData, 2, function(x) x + rnorm(length(x), + sd = abs(x) * refPars$best[12])) + +# Selecting parameters +parSel = c(1:6, 12) + + +## Builidng the likelihood function +likelihood <- function(par, sum = TRUE){ + x = refPars$best + x[parSel] = par + predicted <- VSEM(x[1:11], PAR) + diff = c(predicted[,1:3] - obs[,1:3]) + llValues = dnorm(diff, sd = max(abs(c(predicted[,1:3])),0.0001) * x[12], log = TRUE) + if (sum == False) return(llValues) + else return(sum(llValues)) +} + +# Prior +prior <- createUniformPrior(lower = refPars$lower[parSel], upper = refPars$upper[parSel]) + + +#### +## Definition of the packages and objects that are exported to the cluster. +# These are the objects that are used in the likelihood function. +opts <- list(packages = list("BayesianTools"), variables = list("refPars", "obs", "PAR" ), + dlls = NULL) + + +# Create Bayesian Setup +BSVSEM <- createBayesianSetup(likelihood, prior, best = refPars$best[parSel], + names = rownames(refPars)[parSel], parallel = 2, + parallelOptions = opts)#> Error in get(name, envir = envir): object 'refPars' not found+## The bayesianSetup can now be used in the runMCMC function. +# Note that not all samplers can make use of parallel +# computing. + +# Remove the Bayesian Setup and close the cluster +stopParallel(BSVSEM) +rm(BSVSEM)#> Warning: object 'BSVSEM' not found
Convenience function to create a beta prior
+ + +createBetaPrior(a, b, lower = 0, upper = 1)+ +
a | +shape1 of the beta distribution |
+
---|---|
b | +shape2 of the beta distribution |
+
lower | +lower values for the parameters |
+
upper | +upper values for the parameters |
+
This creates a beta prior, assuming that lower / upper values for parameters are are fixed. The beta is the calculated relative to this lower / upper space.
+ +for details see createPrior
createPriorDensity
+ createPrior
+ createTruncatedNormalPrior
+ createUniformPrior
+ createBayesianSetup
A function for use with plotHist. Creates a matrix representing breaks of a histogram. The matrix will contain the upper bounds, the lower bounds and the frequencies of the breaks in the columns, and the individual breaks in the rows.
+ + +createBreakMat(x, breaks = 15, scale = FALSE)+ +
x | +vector of values |
+
---|---|
breaks | +number of breaks |
+
scale | +logical, if TRUE, the area within the rectangle will be scaled to one (density) |
+
Creates a standardized likelihood class#'
+ + +createLikelihood(likelihood, names = NULL, parallel = F, + catchDuplicates = T, sampler = NULL, parallelOptions = NULL)+ +
likelihood | +Log likelihood density |
+
---|---|
names | +Parameter names (optional) |
+
parallel | +parallelization , either i) no parallelization --> F, ii) native R parallelization --> T / "auto" will select n-1 of your available cores, or provide a number for how many cores to use, or iii) external parallelization --> "external". External means that the likelihood is already able to execute parallel runs in form of a matrix with |
+
catchDuplicates | +Logical, determines whether unique parameter combinations should only be evaluated once. Only used when the likelihood accepts a matrix with parameter as columns. |
+
sampler | +sampler |
+
parallelOptions | +list containing two lists. First "packages" determines the R packages necessary to run the likelihood function. Second "objects" the objects in the global envirnment needed to run the likelihood function (for details see |
+
Convenience function to create an object of class mcmcSamplerList from a list of mcmc samplers
+ + +createMcmcSamplerList(mcmcList)+ +
mcmcList | +a list with each object being an mcmcSampler |
+
---|
Object of class "mcmcSamplerList"
+ + +Creates a standardized posterior class
+ + +createPosterior(prior, likelihood)+ +
prior | +prior class |
+
---|---|
likelihood | +Log likelihood density |
+
Function is internally used in createBayesianSetup
to create a standarized posterior class.
Creates a standardized prior class
+ + +createPrior(density = NULL, sampler = NULL, lower = NULL, upper = NULL, + best = NULL)+ +
density | +Prior density |
+
---|---|
sampler | +Sampling function for density (optional) |
+
lower | +vector with lower bounds of parameters |
+
upper | +vector with upper bounds of parameter |
+
best | +vector with "best" parameter values |
+
This is the general prior generator. It is highly recommended to not only implement the density, but also the sampler function. If this is not done, the user will have to provide explicit starting values for many of the MCMC samplers. Note the existing, more specialized prior function. If your prior can be created by those, they are preferred. Note also that priors can be created from an existing MCMC output from BT, or another MCMC sample, via createPriorDensity
.
min and max truncate, but not re-normalize the prior density (so, if a pdf that integrated to one is truncated, the integral will in general be smaller than one). For MCMC sampling, this doesn't make a difference, but if absolute values of the prior density are a concern, one should provide a truncated density function for the prior.
+ +createPriorDensity
+ createBetaPrior
+ createUniformPrior
+ createTruncatedNormalPrior
+ createBayesianSetup
+# Create a general prior distribution by specifying an arbitrary density function and a +# corresponding sampling function +density = function(par){ + d1 = dunif(par[1], -2,6, log =TRUE) + d2 = dnorm(par[2], mean= 2, sd = 3, log =TRUE) + return(d1 + d2) +} + +# The sampling is optional but recommended because the MCMCs can generate automatic starting +# conditions if this is provided +sampler = function(n=1){ + d1 = runif(n, -2,6) + d2 = rnorm(n, mean= 2, sd = 3) + return(cbind(d1,d2)) +} + +prior <- createPrior(density = density, sampler = sampler, + lower = c(-3,-3), upper = c(3,3), best = NULL) + + +# Use this prior in an MCMC + +ll <- function(x) sum(dnorm(x, log = TRUE)) # multivariate normal ll +bayesianSetup <- createBayesianSetup(likelihood = ll, prior = prior) + +settings = list(iterations = 1000) +out <- runMCMC(bayesianSetup = bayesianSetup, settings = settings)#> Running DEzs-MCMC, chain 1 iteration 300 of 1002 . Current logp -7.903088 -7.501712 -6.397311 . Please wait! Running DEzs-MCMC, chain 1 iteration 600 of 1002 . Current logp -7.536071 -6.169692 -6.187793 . Please wait! Running DEzs-MCMC, chain 1 iteration 900 of 1002 . Current logp -6.154527 -7.003082 -8.155503 . Please wait! Running DEzs-MCMC, chain 1 iteration 1002 of 1002 . Current logp -6.404899 -7.036155 -6.217254 . Please wait!#>+# see ?createPriorDensity for how to create a new prior from this output + +
Fits a density function to a multivariate sample
+ + +createPriorDensity(sampler, method = "multivariate", eps = 1e-10, + lower = NULL, upper = NULL, best = NULL, ...)+ +
sampler | +an object of class BayesianOutput or a matrix |
+
---|---|
method | +method to generate prior - default and currently only option is multivariate |
+
eps | +numerical precision to avoid singularity |
+
lower | +vector with lower bounds of parameter for the new prior, independent of the input sample |
+
upper | +vector with upper bounds of parameter for the new prior, independent of the input sample |
+
best | +vector with "best" values of parameter for the new prior, independent of the input sample |
+
... | +parameters to pass on to the getSample function |
+
createPrior
+ createBetaPrior
+ createTruncatedNormalPrior
+ createUniformPrior
+ createBayesianSetup
+# Create a BayesianSetup +ll <- generateTestDensityMultiNormal(sigma = "no correlation") +bayesianSetup = createBayesianSetup(likelihood = ll, + lower = rep(-10, 3), + upper = rep(10, 3)) + +settings = list(iterations = 2500) +out <- runMCMC(bayesianSetup = bayesianSetup, settings = settings)#> Running DEzs-MCMC, chain 1 iteration 300 of 2502 . Current logp -13.13027 -13.9082 -12.99289 . Please wait! Running DEzs-MCMC, chain 1 iteration 600 of 2502 . Current logp -17.17077 -15.83237 -14.44773 . Please wait! Running DEzs-MCMC, chain 1 iteration 900 of 2502 . Current logp -14.78237 -12.37018 -13.75772 . Please wait! Running DEzs-MCMC, chain 1 iteration 1200 of 2502 . Current logp -13.2786 -13.84704 -14.69235 . Please wait! Running DEzs-MCMC, chain 1 iteration 1500 of 2502 . Current logp -12.24353 -12.61178 -12.44291 . Please wait! Running DEzs-MCMC, chain 1 iteration 1800 of 2502 . Current logp -11.88731 -13.33655 -14.25091 . Please wait! Running DEzs-MCMC, chain 1 iteration 2100 of 2502 . Current logp -13.29007 -13.83644 -13.27353 . Please wait! Running DEzs-MCMC, chain 1 iteration 2400 of 2502 . Current logp -12.2154 -13.3844 -18.34041 . Please wait! Running DEzs-MCMC, chain 1 iteration 2502 of 2502 . Current logp -18.18653 -14.181 -13.19644 . Please wait!#>+ +newPrior = createPriorDensity(out, method = "multivariate", + eps = 1e-10, lower = rep(-10, 3), + upper = rep(10, 3), best = NULL) + +bayesianSetup <- createBayesianSetup(likelihood = ll, prior = newPrior) + +settings = list(iterations = 1000) +out <- runMCMC(bayesianSetup = bayesianSetup, settings = settings)#> Running DEzs-MCMC, chain 1 iteration 300 of 1002 . Current logp -7.545238 -7.446456 -7.128085 . Please wait! Running DEzs-MCMC, chain 1 iteration 600 of 1002 . Current logp -8.492681 -7.42916 -7.616632 . Please wait! Running DEzs-MCMC, chain 1 iteration 900 of 1002 . Current logp -8.09548 -7.647625 -7.012468 . Please wait! Running DEzs-MCMC, chain 1 iteration 1002 of 1002 . Current logp -7.20372 -9.717723 -10.57818 . Please wait!#>+
Factory that creates a proposal generator
+ + +createProposalGenerator(covariance, gibbsProbabilities = NULL, + gibbsWeights = NULL, otherDistribution = NULL, + otherDistributionLocation = NULL, otherDistributionScaled = F, + message = F, method = "chol", scalingFactor = 2.38)+ +
covariance | +covariance matrix. Can also be vector of the sqrt of diagonal elements --> standard deviation |
+
---|---|
gibbsProbabilities | +optional probabilities for the number of parameters to vary in a Metropolis within gibbs style - for 4 parameters, c(1,1,0.5,0) means that at most 3 parameters will be varied, and it is double as likely to vary one or two than varying 3 |
+
gibbsWeights | +optional probabilities for parameters to be varied in a Metropolis within gibbs style - default ist equal weight for all parameters - for 4 parameters, c(1,1,1,100) would mean that if 2 parameters would be selected, parameter 4 would be 100 times more likely to be picked than the others. If 4 is selected, the remaining parameters have equal probability. |
+
otherDistribution | +optional additinal distribution to be mixed with the default multivariate normal. The distribution needs to accept a parameter vector (to allow for the option of making the distribution dependend on the parameter values), but it is still assumed that the change from the current values is returned, not the new absolute values. |
+
otherDistributionLocation | +a vector with 0 and 1, denoting which parameters are modified by the otherDistribution |
+
otherDistributionScaled | +should the other distribution be scaled if gibbs updates are calculated? |
+
message | +print out parameter settings |
+
method | +method for covariance decomposition |
+
scalingFactor | +scaling factor for the proposals |
+
+testMatrix = matrix(rep(c(0,0,0,0), 1000), ncol = 4) +testVector = c(0,0,0,0) + + +##Standard multivariate normal proposal generator + +testGenerator <- createProposalGenerator(covariance = c(1,1,1,1), message = TRUE)#> Proposalgenerator createdcovariance set to: +#> [,1] [,2] [,3] [,4] +#> [1,] 1 0 0 0 +#> [2,] 0 1 0 0 +#> [3,] 0 0 1 0 +#> [4,] 0 0 0 1 +#> covarianceDecomp set to: +#> [,1] [,2] [,3] [,4] +#> [1,] 1 0 0 0 +#> [2,] 0 1 0 0 +#> [3,] 0 0 1 0 +#> [4,] 0 0 0 1 +#> gibbsProbabilities set to: +#> NULL +#> gibbsWeights set to: +#> NULL +#> otherDistribution set to: +#> NULL +#> otherDistributionLocation set to: +#> NULL+methods(class = "proposalGenerator")#> [1] print +#> see '?methods' for accessing help and source codeprint(testGenerator)#> covariance set to: +#> [,1] [,2] [,3] [,4] +#> [1,] 1 0 0 0 +#> [2,] 0 1 0 0 +#> [3,] 0 0 1 0 +#> [4,] 0 0 0 1 +#> covarianceDecomp set to: +#> [,1] [,2] [,3] [,4] +#> [1,] 1 0 0 0 +#> [2,] 0 1 0 0 +#> [3,] 0 0 1 0 +#> [4,] 0 0 0 1 +#> gibbsProbabilities set to: +#> NULL +#> gibbsWeights set to: +#> NULL +#> otherDistribution set to: +#> NULL +#> otherDistributionLocation set to: +#> NULL+x = testGenerator$returnProposal(testVector) +x#> [1] -0.08528815 0.96024687 -0.21403684 0.86855229+x <- testGenerator$returnProposalMatrix(testMatrix) +boxplot(x)+##Changing the covariance +testGenerator$covariance = diag(rep(100,4)) +testGenerator <- testGenerator$updateProposalGenerator(testGenerator, message = TRUE)#> Proposalgenerator settings changedcovariance set to: +#> [,1] [,2] [,3] [,4] +#> [1,] 100 0 0 0 +#> [2,] 0 100 0 0 +#> [3,] 0 0 100 0 +#> [4,] 0 0 0 100 +#> covarianceDecomp set to: +#> [,1] [,2] [,3] [,4] +#> [1,] 10 0 0 0 +#> [2,] 0 10 0 0 +#> [3,] 0 0 10 0 +#> [4,] 0 0 0 10 +#> gibbsProbabilities set to: +#> NULL +#> gibbsWeights set to: +#> NULL +#> otherDistribution set to: +#> NULL +#> otherDistributionLocation set to: +#> NULL+testGenerator$returnProposal(testVector)#> [1] -6.237000 2.271769 0.885039 -3.262444x <- testGenerator$returnProposalMatrix(testMatrix) +boxplot(x)+ +##-Changing the gibbs probabilities / probability to modify 1-n parameters + +testGenerator$gibbsProbabilities = c(1,1,0,0) +testGenerator <- testGenerator$updateProposalGenerator(testGenerator) + +testGenerator$returnProposal(testVector)#> [1] -27.17082 0.00000 -26.53460 0.00000x <- testGenerator$returnProposalMatrix(testMatrix) +boxplot(x)+ +##-Changing the gibbs weights / probability to pick each parameter + +testGenerator$gibbsWeights = c(0.3,0.3,0.3,100) +testGenerator <- testGenerator$updateProposalGenerator(testGenerator) + +testGenerator$returnProposal(testVector)#> [1] 11.78213 0.00000 0.00000 17.64915x <- testGenerator$returnProposalMatrix(testMatrix) +boxplot(x)+ +##-Adding another function + +otherFunction <- function(x) sample.int(10,1) + +testGenerator <- createProposalGenerator( + covariance = c(1,1,1), + otherDistribution = otherFunction, + otherDistributionLocation = c(0,0,0,1), + otherDistributionScaled = TRUE +) + +testGenerator$returnProposal(testVector)#> [1] 0.7176416 1.8980108 -2.8475711 3.5700000x <- testGenerator$returnProposalMatrix(testMatrix) +boxplot(x)table(x[,4])#> +#> 1.19 2.38 3.57 4.76 5.95 7.14 8.33 9.52 10.71 11.9 +#> 129 92 113 97 96 87 89 99 98 100
Convenience function to create an object of class SMCSamplerList from a list of mcmc samplers
+ + +createSmcSamplerList(...)+ +
... | +a list of MCMC samplers |
+
---|
a list of class smcSamplerList with each object being an smcSampler
+ + +Convenience function to create a truncated normal prior
+ + +createTruncatedNormalPrior(mean, sd, lower, upper)+ +
mean | +best estimate for each parameter |
+
---|---|
sd | +sdandard deviation |
+
lower | +vector of lower prior range for all parameters |
+
upper | +vector of upper prior range for all parameters |
+
for details see createPrior
createPriorDensity
+ createPrior
+ createBetaPrior
+ createUniformPrior
+ createBayesianSetup
+prior <- createTruncatedNormalPrior(c(0,0),c(0.4,5), lower = c(-2,-2), upper = c(1,1)) +prior$density(c(2,3))#> [1] -Infprior$density(c(0.2,0.9))#> [1] -1.216469prior$sampler()#> [1] 0.2413041 -1.7412878
Convenience function to create a simple uniform prior distribution
+ + +createUniformPrior(lower, upper, best = NULL)+ +
lower | +vector of lower prior range for all parameters |
+
---|---|
upper | +vector of upper prior range for all parameters |
+
best | +vector with "best" values for all parameters |
+
for details see createPrior
createPriorDensity
, createPrior
, createBetaPrior
, createTruncatedNormalPrior
, createBayesianSetup
+set.seed(1) + +prior <- createUniformPrior(lower = c(0,0), upper = c(0.4,5)) + +# c(2, 3) outside of limits +prior$density(c(2, 3))#> [1] -Inf# -Inf + +# c(0.2, 2) within limits +prior$density(c(0.2, 2))#> [1] -0.6931472# -0.6931472 + + +# sample from prior +prior$sampler()#> [1] 0.2291413 4.5410389# [1] 0.2291413 4.5410389 + + +## the prior object can be passed to createBayesianSetup() + +# log-likelihood density function (needed for createBayesianSetup) +ll <- function(x) sum(dnorm(x, log = TRUE)) + +setup <- createBayesianSetup(prior = prior, likelihood = ll)
Runs Gelman Diagnotics over an BayesianOutput
+ + +gelmanDiagnostics(sampler, thin = "auto", plot = F, ...)+ +
sampler | +an object of class mcmcSampler or mcmcSamplerList |
+
---|---|
thin | +parameter determining the thinning intervall. Either an integer or "auto" (default) for automatic thinning. |
+
plot | +should a Gelman plot be generated |
+
... | +further arguments passed to |
+
The function calls the coda package to calculate Gelman diagnostics and plots
+The original idea is that this function is applied to the outcome of several independent MCMC runs. Technically and practically, it can also be applied to a single MCMC run that has several internal chains, such as DE, DEzs, DREAM, DREAMzs or T-Walk. As argued in ter Braak et al. (2008), the internal chains should be independent after burn-in. While this is likely correct, it also means that they are not completely independent before, and we observed this behavior in the use of the algorithms (i.e. that internal DEzs chains are more similar to each other than the chains of independent DEzs algorithms). A concern is that this non-independence could lead to a failure to detect that the sampler hasn't converged yet. We would therefore recommend to run several DEzs and check convergence with those, instead of running only one.
+ter Braak, Cajo JF, and Jasper A. Vrugt. "Differential evolution Markov chain with snooker updater and fewer chains." Statistics and Computing 18.4 (2008): 435-446.
+ + +Factory to generate a parallel executer of an existing function
+ + +generateParallelExecuter(fun, parallel = F, parallelOptions = list(variables + = "all", packages = "all", dlls = NULL))+ +
fun | +function to be changed to parallel exectution |
+
---|---|
parallel | +should a parallel R cluster be used or not. If set to T, cores will be detected automatically and n-1 of the available n cores of the machine will be used. Alternatively, you can set the number of cores used by hand |
+
parallelOptions | +list containing three lists. First "packages" determines the R packages necessary to run the likelihood function. Second "variables" the objects in the global envirnment needed to run the likelihood function and third "dlls" the DLLs needed to run the likelihood function (see Details). |
+
For parallelization, option T means that an automatic parallelization via R is attempted, or "external", in which case it is assumed that the likelihood is already parallelized. In this case it needs to accept a matrix with parameters as columns. +Further you can specify the packages, objects and DLLs that are exported to the cluster. +By default a copy of your workspace is exported. However, depending on your workspace this can be very inefficient.
+Alternatively you can specify the environments and packages in the likelihood function (e.g. BayesianTools::VSEM() instead of VSEM()).
+ +Can also be used to make functions compatible with library sensitivity
+ + +++testDensityMultiNormal <- generateTestDensityMultiNormal() + + +parDen <- generateParallelExecuter(testDensityMultiNormal)$parallelFun +x = matrix(runif(9,0,1), nrow = 3) +parDen(x)#> [1] -734.8756 -2986.1006 -1098.0516+
Generates a 3 dimensional multivariate normal likelihood function.
+ + +generateTestDensityMultiNormal(mean = c(0, 0, 0), + sigma = "strongcorrelation", sample = F, n = 1, throwErrors = -1)+ +
mean | +vector with the three mean values of the distribution |
+
---|---|
sigma | +either a correlation matrix, or "strongcorrelation", or "no correlation" |
+
sample | +should the function create samples |
+
n | +number of samples to create |
+
throwErrors | +parameter for test purpose. Between 0 and 1 for proportion of errors |
+
3-d multivariate normal density function with mean 2,4,0 and either strong correlation (default), or no correlation.
+ +testDensityBanana
+ testLinearModel
+# sampling from the test function +x = generateTestDensityMultiNormal(sample = TRUE, n = 1000)(1000) +correlationPlot(x)marginalPlot(x)#> Warning: Parameter 'mat' is not of class 'bayesianOutput', set plotPrior to FALSE.+# generating the the density +density = generateTestDensityMultiNormal(sample = FALSE) +density(x[1,])#> [1] 0.2472597
Calculate confidence region from an MCMC or similar sample
+ + +getCredibleIntervals(sampleMatrix, quantiles = c(0.025, 0.975))+ +
sampleMatrix | +matrix of outcomes. Could be parameters or predictions |
+
---|---|
quantiles | +quantiles to be calculated |
+
Creates a DHARMa object
+ + +getDharmaResiduals(model, parMatrix, numSamples, observed, error, plot = TRUE)+ +
model | +function that calculates model predictions for a given parameter vector |
+
---|---|
parMatrix | +a parameter matrix from which the simulations will be generated |
+
numSamples | +the number of samples |
+
observed | +a vector of observed values |
+
error | +function with signature f(mean, par) that generates error expectations from mean model predictions. Par is a vector from the matrix with the parameter samples (full length). f needs to know which of these parameters are parameters of the error function |
+
plot | +logical, determining whether the simulated residuals should be plotted |
+
Returns Metropolis default settings
+ + +getMetropolisDefaultSettings()
+
+
+ Returns possible sampler types
+ + +getPossibleSamplerTypes()
+
+
+ Calculates predictive distribution based on the parameters
+ + +getPredictiveDistribution(parMatrix, model, numSamples = 1000)+ +
parMatrix | +matrix of parameter values |
+
---|---|
model | +model / function to calculate predictions. Outcome should be a vector |
+
numSamples | +number of samples to be drawn |
+
If numSamples is greater than the number of rows in parMatrix, or NULL, or FALSE, or less than 1 all samples in parMatrix will be used.
+ +Calculates Bayesian credible (confidence) and predictive intervals based on parameter sample
+ + +getPredictiveIntervals(parMatrix, model, numSamples = 1000, + quantiles = c(0.025, 0.975), error = NULL)+ +
parMatrix | +matrix of parameter values |
+
---|---|
model | +model / function to calculate predictions. Outcome should be a vector |
+
numSamples | +number of samples to be drawn |
+
quantiles | +quantiles to calculate |
+
error | +function with signature f(mean, par) that generates error expectations from mean model predictions. Par is a vector from the matrix with the parameter samples (full length). f needs to know which of these parameters are parameters of the error function. If supplied, will calculate also predictive intervals additional to credible intervals |
+
If numSamples is greater than the number of rows in parMatrix, or NULL, or FALSE, or less than 1 all samples in parMatrix will be used.
+ +Extracts the sample from a bayesianOutput
+ + +getSample(sampler, parametersOnly = T, coda = F, start = 1, end = NULL, + thin = 1, numSamples = NULL, whichParameters = NULL, + includesProbabilities = F, reportDiagnostics = FALSE, ...)+ +
sampler | +an object of class mcmcSampler, mcmcSamplerList, smcSampler, smcSamplerList, mcmc, mcmc.list, double, numeric |
+
---|---|
parametersOnly | +if F, likelihood, posterior and prior values are also provided in the output |
+
coda | +works only for mcmc classes - provides output as a coda object. Note: if mcmcSamplerList contains mcmc samplers such as DE that have several chains, the internal chains will be collapsed. This may not be the desired behavior for all applications. |
+
start | +for mcmc samplers start value in the chain. For SMC samplers, start particle |
+
end | +for mcmc samplers end value in the chain. For SMC samplers, end particle |
+
thin | +thinning parameter. Either an integer determining the thinning intervall (default is 1) or "auto" for automatic thinning. |
+
numSamples | +sample size (only used if thin = 1). If you want to use numSamples set thin to 1. |
+
whichParameters | +possibility to select parameters by index |
+
includesProbabilities | +applies only to getSample.Matrix. logical, determining whether probabilities should be included in the result. |
+
reportDiagnostics | +logical, determines whether settings should be included in the output |
+
... | +further arguments |
+
If thin is greater than the total number of samples in the sampler object the first and the last element (of each chain if a sampler with multiples chains is used) are sampled. If numSamples is greater than the total number of samples all samples are selected. In both cases a warning is displayed.
+If thin and numSamples is passed, the function will use the thin argument if it is valid and greater than 1, else numSamples will be used.
+ + ++ll = function(x) sum(dnorm(x, log = TRUE)) + +setup = createBayesianSetup(ll, lower = c(-10,-10), upper = c(10,10)) + +settings = list(nrChains = 2, iterations = 1000) +out <- runMCMC(bayesianSetup = setup, sampler = "DEzs", settings = settings)#> Running DEzs-MCMC, chain 1 iteration 300 of 1002 . Current logp -14.26114 -8.366528 -8.935835 . Please wait! Running DEzs-MCMC, chain 1 iteration 600 of 1002 . Current logp -9.937957 -8.081307 -7.929941 . Please wait! Running DEzs-MCMC, chain 1 iteration 900 of 1002 . Current logp -8.117279 -8.05287 -9.947761 . Please wait! Running DEzs-MCMC, chain 1 iteration 1002 of 1002 . Current logp -8.392833 -10.03726 -9.426219 . Please wait!#>#> Running DEzs-MCMC, chain 2 iteration 300 of 1002 . Current logp -7.991262 -8.611373 -7.993276 . Please wait! Running DEzs-MCMC, chain 2 iteration 600 of 1002 . Current logp -8.018305 -13.25595 -8.069903 . Please wait! Running DEzs-MCMC, chain 2 iteration 900 of 1002 . Current logp -8.563249 -7.926249 -8.226997 . Please wait! Running DEzs-MCMC, chain 2 iteration 1002 of 1002 . Current logp -7.835189 -8.293727 -8.162822 . Please wait!#>+# population MCMCs divide the interations by the number of internal chains, +# so the end of the 3 chains is 1000/3 = 334 +sample <- getSample(out, start = 100, end = 334, thin = 10) + +# sampling with number of samples instead of thinning and +# returning a coda object +sample <- getSample(out, start = 100, numSamples = 60, coda = TRUE) +plot(sample)+ +# MCMC with a single chain: +settings_2 <- list(nrChains = 1, iterations = 1000) +out_2 <- runMCMC(setup, sampler = "Metropolis", settings = settings_2)#> BT runMCMC: trying to find optimal start and covariance values#>#> Running Metropolis-MCMC, chain 1 iteration 100 of 1000 . Current logp: -9.532038 Please wait! Running Metropolis-MCMC, chain 1 iteration 200 of 1000 . Current logp: -8.668119 Please wait! Running Metropolis-MCMC, chain 1 iteration 300 of 1000 . Current logp: -7.900279 Please wait! Running Metropolis-MCMC, chain 1 iteration 400 of 1000 . Current logp: -8.629034 Please wait! Running Metropolis-MCMC, chain 1 iteration 500 of 1000 . Current logp: -8.510623 Please wait! Running Metropolis-MCMC, chain 1 iteration 600 of 1000 . Current logp: -9.273203 Please wait! Running Metropolis-MCMC, chain 1 iteration 700 of 1000 . Current logp: -7.956503 Please wait! Running Metropolis-MCMC, chain 1 iteration 800 of 1000 . Current logp: -8.932149 Please wait! Running Metropolis-MCMC, chain 1 iteration 900 of 1000 . Current logp: -8.071972 Please wait! Running Metropolis-MCMC, chain 1 iteration 1000 of 1000 . Current logp: -8.17257 Please wait!#>sample_2 <- getSample(out_2, numSamples = 100)
Calculate posterior volume
+ + +getVolume(sampler, prior = F, method = "MVN", ...)+ +
sampler | +an object of superclass bayesianOutput or any other class that has the getSample function implemented (e.g. Matrix) |
+
---|---|
prior | +schould also prior volume be calculated |
+
method | +method for volume estimation. Currently, the only option is "MVN" |
+
... | +additional parameters to pass on to the |
+
The idea of this function is to provide an estimate of the "posterior volume", i.e. how "broad" the posterior is. One potential application is to the overall reduction of parametric uncertainty between different data types, or between prior and posterior.
+Implemented methods for volume estimation:
+Option "MVN" - in this option, the volume is calculated as the determinant of the covariance matrix of the prior / posterior sample.
+ + ++bayesianSetup = createBayesianSetup( + likelihood = generateTestDensityMultiNormal(sigma = "no correlation"), + lower = rep(-10, 3), upper = rep(10, 3)) + +out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", + settings = list(iterations = 2000, message = FALSE)) + +getVolume(out, prior = TRUE)#> $priorVol +#> [1] 38943.63 +#> +#> $postVol +#> [1] 1.035675 +#>+bayesianSetup = createBayesianSetup( + likelihood = generateTestDensityMultiNormal(sigma = "strongcorrelation"), + lower = rep(-10, 3), upper = rep(10, 3)) + +out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", + settings = list(iterations = 2000, message = FALSE)) + +getVolume(out, prior = TRUE)#> $priorVol +#> [1] 38388.73 +#> +#> $postVol +#> [1] 2.82539e-06 +#>+
Internal function to plot histograms for the marginalPlot function
+ + +histMarginal(x, at = 0, .range = NULL, breaks = 15, col = c("#FF5000C0", + "#4682B4A0"), border = c("black", "black"), main = "", dens = FALSE, + densCol = c("black", "black"), densLwd = c(2, 2), densLty = c(2, 2), + add = TRUE, res = 500)+ +
x | +list of vectors. Each vector will be plotted as a separate histogram. |
+
---|---|
at | +y position of the histograms |
+
.range | +maximum height of the histogram. If NULL, will be determined automatically |
+
breaks | +integer, determining the number of breaks |
+
col | +vector of colors determining histogram colors |
+
border | +vector of colors determining histogram border colors |
+
main | +Character, determining the title of the plot |
+
dens | +Logical, determining whether a density plot should be plotted additionally to the histogram |
+
densCol | +vector of colors for the density plots |
+
densLwd | +line width of the density plot |
+
densLty | +line type of the density plot |
+
add | +Logical, determining whether the histogram should be added to and existing plot window. |
+
res | +resolution of the density overlay |
+
AR1 type likelihood function
+ + +likelihoodAR1(predicted, observed, sd, a)+ +
predicted | +vector of predicted values |
+
---|---|
observed | +vector of observed values |
+
sd | +standard deviation of the iid normal likelihood |
+
a | +temporal correlation in the AR1 model |
+
The AR1 model considers the process: y(t) = a y(t-1) + E e = i.i.d. N(0,sd) |a| < 1 At the moment, no NAs are allowed in the time series.
+ + +Normal / Gaussian Likelihood function
+ + +likelihoodIidNormal(predicted, observed, sd)+ +
predicted | +vector of predicted values |
+
---|---|
observed | +vector of observed values |
+
sd | +standard deviation of the i.i.d. normal likelihood |
+
Calcluated the marginal likelihood from a set of MCMC samples
+ + +marginalLikelihood(sampler, numSamples = 1000, method = "Chib", ...)+ +
sampler | +an object that implements the getSample function, i.e. a mcmc / smc Sampler (list) |
+
---|---|
numSamples | +number of samples to use. How this works, and if it requires recalculating the likelihood, depends on the method |
+
method | +method to choose. Currently available are "Chib" (default), the harmonic mean "HM", sampling from the prior "prior", and bridge sampling "Bridge". See details |
+
... | +further arguments passed to |
+
The function currently implements four ways to calculate the marginal likelihood. + The recommended way is the method "Chib" (Chib and Jeliazkov, 2001). which is based on MCMC samples, but performs additional calculations. + Despite being the current recommendation, note there are some numeric issues with this algorithm that may limit reliability for larger dimensions.
+The harmonic mean approximation, + is implemented only for comparison. Note that the method is numerically + unrealiable and usually should not be used.
+The third method is simply sampling from the prior. While in principle unbiased, + it will only converge for a large number of samples, and is therefore + numerically inefficient.
+The Bridge method uses bridge sampling as implemented in the R package "bridgesampling".
+ +Chib, Siddhartha, and Ivan Jeliazkov. "Marginal likelihood from the Metropolis-Hastings output." Journal of the American Statistical Association 96.453 (2001): 270-281.
+ ++# Harmonic mean works OK for a low-dim case with + +likelihood <- function(x) sum(dnorm(x, log = TRUE)) +prior = createUniformPrior(lower = rep(-1,2), upper = rep(1,2)) +bayesianSetup <- createBayesianSetup(likelihood = likelihood, prior = prior) +out = runMCMC(bayesianSetup = bayesianSetup, settings = list(iterations = 5000))#> Running DEzs-MCMC, chain 1 iteration 300 of 5001 . Current logp -3.231946 -3.22838 -3.702234 . Please wait! Running DEzs-MCMC, chain 1 iteration 600 of 5001 . Current logp -3.650859 -3.465972 -3.512521 . Please wait! Running DEzs-MCMC, chain 1 iteration 900 of 5001 . Current logp -3.609262 -3.471144 -3.339216 . Please wait! Running DEzs-MCMC, chain 1 iteration 1200 of 5001 . Current logp -3.58876 -3.691328 -3.307711 . Please wait! Running DEzs-MCMC, chain 1 iteration 1500 of 5001 . Current logp -3.821757 -3.323898 -3.385563 . Please wait! Running DEzs-MCMC, chain 1 iteration 1800 of 5001 . Current logp -3.23847 -3.725506 -3.288343 . Please wait! Running DEzs-MCMC, chain 1 iteration 2100 of 5001 . Current logp -3.394055 -3.738272 -3.327836 . Please wait! Running DEzs-MCMC, chain 1 iteration 2400 of 5001 . Current logp -3.548363 -3.569064 -3.510509 . Please wait! Running DEzs-MCMC, chain 1 iteration 2700 of 5001 . Current logp -3.759792 -3.226294 -3.723042 . Please wait! Running DEzs-MCMC, chain 1 iteration 3000 of 5001 . Current logp -3.645238 -3.313713 -3.421209 . Please wait! Running DEzs-MCMC, chain 1 iteration 3300 of 5001 . Current logp -3.551052 -3.465299 -3.66877 . Please wait! Running DEzs-MCMC, chain 1 iteration 3600 of 5001 . Current logp -3.272982 -3.333179 -3.631589 . Please wait! Running DEzs-MCMC, chain 1 iteration 3900 of 5001 . Current logp -3.250819 -3.700023 -3.709891 . Please wait! Running DEzs-MCMC, chain 1 iteration 4200 of 5001 . Current logp -3.312153 -3.270711 -3.244133 . Please wait! Running DEzs-MCMC, chain 1 iteration 4500 of 5001 . Current logp -3.728044 -3.569831 -4.068991 . Please wait! Running DEzs-MCMC, chain 1 iteration 4800 of 5001 . Current logp -3.370315 -3.700947 -3.653093 . Please wait! Running DEzs-MCMC, chain 1 iteration 5001 of 5001 . Current logp -3.400083 -3.660208 -3.874071 . Please wait!#>+plot(out)+ +marginalLikelihood(out, numSamples = 500)[[1]]#> Warning: Note to the user: be aware that marginal likelihood calculations are notoriously prone to numerical stability issues. Especially in high-dimensional parameter spaces, there is no guarantee that the algorithms implemented in this function converge in all cases. Proceed at your own risk!#> [1] -2.045501marginalLikelihood(out, method = "HM", numSamples = 500)[[1]]#> Warning: Note to the user: be aware that marginal likelihood calculations are notoriously prone to numerical stability issues. Especially in high-dimensional parameter spaces, there is no guarantee that the algorithms implemented in this function converge in all cases. Proceed at your own risk!#> [1] -2.148932marginalLikelihood(out, method = "Prior", numSamples = 500)[[1]]#> Warning: Note to the user: be aware that marginal likelihood calculations are notoriously prone to numerical stability issues. Especially in high-dimensional parameter spaces, there is no guarantee that the algorithms implemented in this function converge in all cases. Proceed at your own risk!#> [1] -2.144425marginalLikelihood(out, method = "Bridge", numSamples = 500)[[1]]#> Iteration: 1 +#> Iteration: 2 +#> Iteration: 3 +#> Iteration: 4 +#> Iteration: 5 +#> Iteration: 6 +#> Iteration: 7#> Warning: Note to the user: be aware that marginal likelihood calculations are notoriously prone to numerical stability issues. Especially in high-dimensional parameter spaces, there is no guarantee that the algorithms implemented in this function converge in all cases. Proceed at your own risk!#> [1] -2.108433+# True marginal likelihood (brute force approximation) + +marginalLikelihood(out, method = "Prior", numSamples = 10000)[[1]]#> Warning: Note to the user: be aware that marginal likelihood calculations are notoriously prone to numerical stability issues. Especially in high-dimensional parameter spaces, there is no guarantee that the algorithms implemented in this function converge in all cases. Proceed at your own risk!#> [1] -2.152117+ +# Harmonic mean goes totally wrong for higher dimensions - wide prior. +# Same goes for standard bridge sampling. +# Could also be a problem of numeric stability of the implementation + +likelihood <- function(x) sum(dnorm(x, log = TRUE)) +prior = createUniformPrior(lower = rep(-10,3), upper = rep(10,3)) +bayesianSetup <- createBayesianSetup(likelihood = likelihood, prior = prior) +out = runMCMC(bayesianSetup = bayesianSetup, settings = list(iterations = 5000))#> Running DEzs-MCMC, chain 1 iteration 300 of 5001 . Current logp -12.3117 -16.32338 -13.25312 . Please wait! Running DEzs-MCMC, chain 1 iteration 600 of 5001 . Current logp -13.15733 -12.68131 -14.6254 . Please wait! Running DEzs-MCMC, chain 1 iteration 900 of 5001 . Current logp -12.55761 -12.41566 -12.51362 . Please wait! Running DEzs-MCMC, chain 1 iteration 1200 of 5001 . Current logp -11.8798 -12.53557 -14.97952 . Please wait! Running DEzs-MCMC, chain 1 iteration 1500 of 5001 . Current logp -12.97593 -12.59545 -13.68896 . Please wait! Running DEzs-MCMC, chain 1 iteration 1800 of 5001 . Current logp -12.97854 -13.25449 -12.38997 . Please wait! Running DEzs-MCMC, chain 1 iteration 2100 of 5001 . Current logp -12.13808 -13.44797 -14.34677 . Please wait! Running DEzs-MCMC, chain 1 iteration 2400 of 5001 . Current logp -11.90999 -12.65168 -14.10616 . Please wait! Running DEzs-MCMC, chain 1 iteration 2700 of 5001 . Current logp -13.22604 -13.94732 -12.74675 . Please wait! Running DEzs-MCMC, chain 1 iteration 3000 of 5001 . Current logp -12.94641 -11.99008 -12.02827 . Please wait! Running DEzs-MCMC, chain 1 iteration 3300 of 5001 . Current logp -11.80098 -13.23051 -16.26073 . Please wait! Running DEzs-MCMC, chain 1 iteration 3600 of 5001 . Current logp -12.41389 -12.22136 -14.67628 . Please wait! Running DEzs-MCMC, chain 1 iteration 3900 of 5001 . Current logp -14.31968 -14.49065 -12.45884 . Please wait! Running DEzs-MCMC, chain 1 iteration 4200 of 5001 . Current logp -12.56762 -14.59702 -12.7141 . Please wait! Running DEzs-MCMC, chain 1 iteration 4500 of 5001 . Current logp -12.98071 -13.00499 -13.57345 . Please wait! Running DEzs-MCMC, chain 1 iteration 4800 of 5001 . Current logp -12.01666 -12.73201 -12.8051 . Please wait! Running DEzs-MCMC, chain 1 iteration 5001 of 5001 . Current logp -14.10143 -12.14869 -14.50727 . Please wait!#>+plot(out)+marginalLikelihood(out, numSamples = 500)[[1]]#> Warning: Note to the user: be aware that marginal likelihood calculations are notoriously prone to numerical stability issues. Especially in high-dimensional parameter spaces, there is no guarantee that the algorithms implemented in this function converge in all cases. Proceed at your own risk!#> [1] -8.986955marginalLikelihood(out, method = "HM", numSamples = 500)[[1]]#> Warning: Note to the user: be aware that marginal likelihood calculations are notoriously prone to numerical stability issues. Especially in high-dimensional parameter spaces, there is no guarantee that the algorithms implemented in this function converge in all cases. Proceed at your own risk!#> [1] -38.51889marginalLikelihood(out, method = "Prior", numSamples = 500)[[1]]#> Warning: Note to the user: be aware that marginal likelihood calculations are notoriously prone to numerical stability issues. Especially in high-dimensional parameter spaces, there is no guarantee that the algorithms implemented in this function converge in all cases. Proceed at your own risk!#> [1] -9.173093marginalLikelihood(out, method = "Bridge", numSamples = 500)[[1]]#> Iteration: 1 +#> Iteration: 2 +#> Iteration: 3 +#> Iteration: 4 +#> Iteration: 5 +#> Iteration: 6 +#> Iteration: 7 +#> Iteration: 8 +#> Iteration: 9 +#> Iteration: 10 +#> Iteration: 11 +#> Iteration: 12 +#> Iteration: 13 +#> Iteration: 14 +#> Iteration: 15 +#> Iteration: 16 +#> Iteration: 17 +#> Iteration: 18#> Warning: Note to the user: be aware that marginal likelihood calculations are notoriously prone to numerical stability issues. Especially in high-dimensional parameter spaces, there is no guarantee that the algorithms implemented in this function converge in all cases. Proceed at your own risk!#> [1] -41.42371+# True marginal likelihood (brute force approximation) + +marginalLikelihood(out, method = "Prior", numSamples = 10000)[[1]]#> Warning: Note to the user: be aware that marginal likelihood calculations are notoriously prone to numerical stability issues. Especially in high-dimensional parameter spaces, there is no guarantee that the algorithms implemented in this function converge in all cases. Proceed at your own risk!#> [1] -9.229592+ +
Plot MCMC marginals
+ + +marginalPlot(mat, thin = "auto", scale = NULL, best = NULL, + histogram = TRUE, plotPrior = TRUE, priorTop = FALSE, + nDrawsPrior = 1000, breaks = 15, res = 500, singlePanel = FALSE, + dens = TRUE, col = c("#FF5000D0", "#4682B4A0"), lwd = par("lwd"), ...)+ +
mat | +object of class "bayesianOutput" or a matrix or data frame of variables |
+
---|---|
thin | +thinning of the matrix to make things faster. Default is to thin to 5000 |
+
scale | +should the results be scaled. Value can be either NULL (no scaling), T, or a matrix with upper / lower bounds as columns. If set to T, attempts to retrieve the scaling from the input object mat (requires that this is of class BayesianOutput) |
+
best | +if provided, will draw points at the given values (to display true / default parameter values). Value can be either NULL (no drawing), a vector with values, or T, in which case the function will attempt to retrieve the values from a BayesianOutput |
+
histogram | +Logical, determining whether a violin plot or a histogram should be plotted |
+
plotPrior | +Logical, determining whether the prior should be plotted in addition to the posteriors. Only applicable if mat is an object of class "bayesianOutput" |
+
priorTop | +Logical, determining whether the prior should be plotted top (TRUE) or bottom (FALSE) |
+
nDrawsPrior | +Integer, number of draws from the prior, when plotPrior is active |
+
breaks | +Integer, number of histogram breaks if histogram is set to TRUE |
+
res | +resolution parameter for violinPlot, determining how many descrete points should be used for the density kernel. |
+
singlePanel | +Logical, determining whether all histograms/violins should be plotted in a single plot panel or in separate panels. |
+
dens | +Logical, determining wheter an density overlay should be plotted when 'histogram' is TRUE |
+
col | +vector of colors for posterior and prior |
+
lwd | +line width of the violin plots |
+
... | +additional parameters to pass on to the |
+
+#> Warning: Parameter 'mat' is not of class 'bayesianOutput', set plotPrior to FALSE.
Run multiple chains
+ + +mcmcMultipleChains(bayesianSetup, settings, sampler)+ +
bayesianSetup | +Object of class "BayesianSetup" |
+
---|---|
settings | +list with settings for sampler |
+
sampler | +character, either "Metropolis" or "DE" |
+
list containing the single runs ($sampler) and the chains in a coda::mcmc.list ($mcmc.list)
+ + +This function is deprecated and will be removed by v0.2.
+ + +createMixWithDefaults(pars, defaults, locations)+ +
pars | +vector with new parameter values |
+
---|---|
defaults | +vector with defaukt parameter values |
+
locations | +indices of the new parameter values |
+
A simple custom histogram plotting function
+ + +plotHist(x, at = NULL, col = "orangered", border = "black")+ +
x | +matrix (e.g. constructed with |
+
---|---|
at | +y position of the histogram. If NULL a new plot will be generated. |
+
col | +color of the histogram |
+
border | +border color |
+
Performs a one-factor-at-a-time sensitivity analysis for the posterior of a given bayesianSetup within the prior range.
+ + +plotSensitivity(bayesianSetup, selection = NULL)+ +
bayesianSetup | +An object of class BayesianSetup |
+
---|---|
selection | +indices of selected parameters |
+
This function can also be used for sensitivity analysis of an arbitrary output - just create a BayesianSetup with this output.
+ + +Plots a time series, with the option to include confidence and prediction band
+ + +plotTimeSeries(observed = NULL, predicted = NULL, x = NULL, + confidenceBand = NULL, predictionBand = NULL, xlab = "Time", + ylab = "Observed / predicted values", ...)+ +
observed | +observed values |
+
---|---|
predicted | +predicted values |
+
x | +optional values for x axis (time) |
+
confidenceBand | +matrix with confidenceBand |
+
predictionBand | +matrix with predictionBand |
+
xlab | +a title for the x axis |
+
ylab | +a title for the y axis |
+
... | +further arguments passed to |
+
plotTimeSeriesResults
+ marginalPlot
+ tracePlot
+ correlationPlot
+# Create time series +ts <- VSEMcreatePAR(1:100) + +# create fake "predictions" +pred <- ts + rnorm(length(ts), mean = 0, sd = 2) + +# plot time series +par(mfrow=c(1,2)) + +plotTimeSeries(observed = ts, main="Observed") +plotTimeSeries(observed = ts, predicted = pred, main = "Observed and predicted")+par(mfrow=c(1,1))
Creates a time series plot typical for an MCMC / SMC fit
+ + +plotTimeSeriesResults(sampler, model, observed, error = NULL, + plotResiduals = TRUE, start = 1, prior = FALSE, ...)+ +
sampler | +Either a) a matrix b) an MCMC object (list or not), or c) an SMC object |
+
---|---|
model | +function that calculates model predictions for a given parameter vector |
+
observed | +observed values |
+
error | +function with signature f(mean, par) that generates error expectations from mean model predictions. Par is a vector from the matrix with the parameter samples (full length). f needs to know which of these parameters are parameters of the error function |
+
plotResiduals | +logical determining whether residuals should be plotted |
+
start | +numeric start value for the plot (see |
+
prior | +if a prior sampler is implemented, setting this parameter to TRUE will draw model parameters from the prior instead of the posterior distribution |
+
... | +further arguments passed to |
+
Rescales values in the interval "from" (lower, upper) to the new interval "to" (lower, upper).
+ + +rescale(x, from, to)+ +
x | +Vector of values |
+
---|---|
from | +vector, interval of which x are elements. from[1] must be the lower, from[2] the upper bound. |
+
to | +vector, interval to which the elements should be scaled. to[1] must be the lower, to[2] the upper bound. |
+
Main wrapper function to start MCMCs, particle MCMCs and SMCs
+ + +runMCMC(bayesianSetup, sampler = "DEzs", settings = NULL)+ +
bayesianSetup | +either one of a) an object of class BayesianSetup with prior and likelihood function (recommended, see |
+
---|---|
sampler | +sampling algorithm to be run. Default is DEzs. Options are "Metropolis", "AM", "DR", "DRAM", "DE", "DEzs", "DREAM", "DREAMzs", "SMC". For details see the help of the individual functions. |
+
settings | +list with settings for each sampler (see help of sampler for details). If a setting is not provided, defaults (see |
+
The function returns an object of class mcmcSampler (if one chain is run) or mcmcSamplerList. Both have the superclass bayesianOutput. It is possible to extract the samples as a coda object or matrix with getSample
.
+It is also possible to summarize the posterior as a new prior via createPriorDensity
.
The runMCMC function can be started with either one of a) an object of class BayesianSetup with prior and likelihood function (recommended, see createBayesianSetup
), b) a log posterior or other target function, or c) an object of class BayesianOutput created by runMCMC. The latter allows to continue a previous MCMC run. If a bayesianSetup is provided, check if appropriate parallization options are used - many samplers can make use of parallelization if this option is activated when the class is created.
For details about the different MCMC samplers, make sure you have read the Vignette (run vignette("BayesianTools", package="BayesianTools"). Also, see Metropolis
for Metropolis based samplers, DE
and DEzs
for standard differential evolution samplers, DREAM
and DREAMzs
for DREAM sampler, Twalk
for the Twalk sampler, and smcSampler
for rejection and Sequential Monte Carlo sampling.
The samplers "AM", "DR", and "DRAM" are special cases of the "Metropolis" sampler and are shortcuts for predefined settings ("AM": adapt=TRUE; "DR": DRlevels=2; "DRAM": adapt=True, DRlevels=2).
+The settings list allows to change the settings for the MCMC samplers and some other options. For the MCMC sampler settings, see their help files. Global options that apply for all MCMC samplers are: iterations (number of MCMC iterations), and nrChains (number of chains to run). Note that running several chains is not done in parallel, so if time is an issue it will be better to run the MCMCs individually and then combine them via createMcmcSamplerList
into one joint object.
Startvalues: all samplers allow to provide explicit startvalues. If startvalues are not provided, they are sampled from the prior. Usually, this is a good choice, so don't feel compelled to provide startvalues.
+Note that DE and DREAM variants as well as SMC and T-walk require a population to start, which should be provided as a matrix. Default (NULL) sets the population size for DE to 3 x dimensions of parameters, for DREAM to 2 x dimensions of parameters and for DEzs and DREAMzs to three, sampled from the prior. Note also that the zs variants of DE and DREAM require two populations, the current population and the z matrix (a kind of memory) - if you want to set both, provide a list with startvalue$X and startvalue$Z.
+Startvalues for sampling with nrChains > 1 : if you want to provide different start values for the different chains, provide them as a list
+ ++## Generate a test likelihood function. +ll <- generateTestDensityMultiNormal(sigma = "no correlation") + +## Create a BayesianSetup object from the likelihood +## is the recommended way of using the runMCMC() function. +bayesianSetup <- createBayesianSetup(likelihood = ll, lower = rep(-10, 3), upper = rep(10, 3)) + +## Finally we can run the sampler and have a look +settings = list(iterations = 1000, adapt = FALSE) +out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings)#> BT runMCMC: trying to find optimal start and covariance values#>#> Running Metropolis-MCMC, chain 1 iteration 100 of 1000 . Current logp: -12.95988 Please wait! Running Metropolis-MCMC, chain 1 iteration 200 of 1000 . Current logp: -13.51121 Please wait! Running Metropolis-MCMC, chain 1 iteration 300 of 1000 . Current logp: -13.21691 Please wait! Running Metropolis-MCMC, chain 1 iteration 400 of 1000 . Current logp: -18.70323 Please wait! Running Metropolis-MCMC, chain 1 iteration 500 of 1000 . Current logp: -12.16547 Please wait! Running Metropolis-MCMC, chain 1 iteration 600 of 1000 . Current logp: -13.34294 Please wait! Running Metropolis-MCMC, chain 1 iteration 700 of 1000 . Current logp: -11.93483 Please wait! Running Metropolis-MCMC, chain 1 iteration 800 of 1000 . Current logp: -12.19979 Please wait! Running Metropolis-MCMC, chain 1 iteration 900 of 1000 . Current logp: -12.35007 Please wait! Running Metropolis-MCMC, chain 1 iteration 1000 of 1000 . Current logp: -13.35738 Please wait!#>+## out is of class bayesianOutput. There are various standard functions +# implemented for this output + +plot(out)correlationPlot(out)marginalPlot(out)summary(out)#> # # # # # # # # # # # # # # # # # # # # # # # # # +#> ## MCMC chain summary ## +#> # # # # # # # # # # # # # # # # # # # # # # # # # +#> +#> # MCMC sampler: Metropolis +#> # Nr. Chains: 1 +#> # Iterations per chain: 1000 +#> # Rejection rate: 0.655 +#> # Effective sample size: 122 +#> # Runtime: 0.73 sec. +#> +#> # Parameters +#> MAP 2.5% median 97.5% +#> par 1 -0.001 -2.038 0.023 2.332 +#> par 2 -0.001 -1.966 0.003 1.948 +#> par 3 0.000 -1.826 0.146 2.202 +#> +#> ## DIC: 12.082 +#> ## Convergence +#> Gelman Rubin multivariate psrf: Only one chain; convergence cannot be determined! +#> +#> ## Correlations +#> par 1 par 2 par 3 +#> par 1 1.000 0.039 0.041 +#> par 2 0.039 1.000 -0.044 +#> par 3 0.041 -0.044 1.000+## additionally, you can return the sample as a coda object, and make use of the coda functions +# for plotting and analysis + +codaObject = getSample(out, start = 500, coda = TRUE)
Gets n equally spaced samples (rows) from a matrix or vector
+ + +sampleEquallySpaced(x, numSamples)+ +
x | +matrix or vector |
+
---|---|
numSamples | +number of samples (rows) to be drawn |
+
Gets n equally spaced samples (rows) from a matrix and returns a new matrix (or vector) containing those samples
+ + +Function to sample with cobinations of the basic Metropolis-Hastings MCMC algorithm (Metropolis et al., 1953), a variation of the adaptive Metropolis MCMC (Haario et al., 2001), the delayed rejection algorithm (Tierney & Mira, 1999), and the delayed rejection adaptive Metropolis algorithm (DRAM, Haario et al), and the Metropolis within Gibbs
+ + +sampleMetropolis(mcmcSampler, iterations)+ +
mcmcSampler | +an mcmcSampler |
+
---|---|
iterations | +iterations |
+
Help function to find starvalues and proposalGenerator settings
+ + +setupStartProposal(proposalGenerator = NULL, bayesianSetup, settings)+ +
proposalGenerator | +proposal generator |
+
---|---|
bayesianSetup | +either an object of class bayesianSetup created by |
+
settings | +list with settings |
+
Sequential Monte Carlo Sampler
+ + +smcSampler(bayesianSetup, initialParticles = 1000, iterations = 10, + resampling = T, resamplingSteps = 2, proposal = NULL, adaptive = T, + proposalScale = 0.5, likelihoodWeights = NULL)+ +
bayesianSetup | +either an object of class bayesianSetup created by |
+
---|---|
initialParticles | +initial particles - either a draw from the prior, provided as a matrix with the single parameters as columns and each row being one particle (parameter vector), or a numeric value with the number of desired particles. In this case, the sampling option must be provided in the prior of the BayesianSetup. |
+
iterations | +number of iterations |
+
resampling | +if new particles should be created at each iteration |
+
resamplingSteps | +how many resampling (MCMC) steps between the iterations |
+
proposal | +optional proposal class |
+
adaptive | +should the covariance of the proposal be adapted during sampling |
+
proposalScale | +scaling factor for the proposal generation. Can be adapted if there is too much / too little rejection |
+
likelihoodWeights | +NULL, numeric vector, or function. The sum of all weights must be one, and the length of the vector must be equal to iterations. The value of the nth element determines the weight for the likelihood in the nth iteration. |
+
The sampler can be used for rejection sampling as well as for sequential Monte Carlo. For the former case set the iterations to one.
+ +The SMC currently assumes that the initial particle is sampled from the prior. If a better initial estimate of the posterior distribution is available, this the sampler should be modified to include this. Currently, however, this is not included in the code, so the appropriate adjustments have to be done by hand.
+ + ++## Example for the use of SMC +# First we need a bayesianSetup - SMC makes most sense if we can for demonstration, +# we'll write a function that puts out the number of model calls + +MultiNomialNoCor <- generateTestDensityMultiNormal(sigma = "no correlation") + +parallelLL <- function(parMatrix){ + print(paste("Calling likelihood with", nrow(parMatrix), "parameter combinations")) + out = apply(parMatrix, 1, MultiNomialNoCor) + return(out) +} + +bayesianSetup <- createBayesianSetup(likelihood = parallelLL, lower = rep(-10, 3), + upper = rep(10, 3), parallel = "external") + +# Defining settings for the sampler +# First we use the sampler for rejection sampling +settings <- list(initialParticles = 1000, iterations = 1, resampling = FALSE) + +# Running the sampler +out1 <- runMCMC(bayesianSetup = bayesianSetup, sampler = "SMC", settings = settings)#> [1] 1 +#> [1] "Calling likelihood with 1000 parameter combinations" +#> [1] 2.984645#> [1] 26#>#plot(out1) + + +# Now for sequential Monte Carlo +settings <- list(initialParticles = 100, iterations = 5, resamplingSteps = 1) +out2 <- runMCMC(bayesianSetup = bayesianSetup, sampler = "SMC", settings = settings)#> [1] 0.2 0.2 0.2 0.2 0.2 +#> [1] "Calling likelihood with 100 parameter combinations" +#> [1] 3.31721#> [1] 17 +#> [1] "Calling likelihood with 2 parameter combinations" +#> [1] "Calling likelihood with 100 parameter combinations" +#> [1] 47.89169#> [1] 49 +#> [1] "Calling likelihood with 2 parameter combinations" +#> [1] "Calling likelihood with 100 parameter combinations" +#> [1] 82.96134#> [1] 66 +#> [1] "Calling likelihood with 3 parameter combinations" +#> [1] "Calling likelihood with 100 parameter combinations" +#> [1] 89.32615#> [1] 63 +#> [1] "Calling likelihood with 100 parameter combinations" +#> [1] 95.49201#> [1] 63 +#> [1] "Calling likelihood with parameter combinations" +#> Parameter values -7.49454366195938 -6.57156059686812 -7.17224688752869#> Warning: Problem encountered in the calculation of the likelihood with parameter -7.49454366195938-6.57156059686812-7.17224688752869 +#> Error message wasError in apply(parMatrix, 1, MultiNomialNoCor): dim(X) must have a positive length +#> +#> set result of the parameter evaluation to -Inf ParaeterValues#>#plot(out2) + +
Banana-shaped density function
+ + +testDensityBanana(p)+ +
p | +2-dim parameter vector |
+
---|
inspired from package FMEmcmc, seems to go back to Laine M (2008). Adaptive MCMC Methods with Applications in Environmental and Models. Finnish Meteorological Institute Contributions 69. ISBN 978-951-697-662-7.
+ +3d Mutivariate Normal likelihood
+ + +testDensityMultiNormal(x, sigma = "strongcorrelation")+ +
x | +a parameter vector of arbitrary length |
+
---|---|
sigma | +either a correlation matrix, or "strongcorrelation", or "no correlation" |
+
Fake model, returns a ax + b linear response to 2-param vector
+ + +testLinearModel(x, env = NULL)+ +
x | +2-dim parameter vector |
+
---|---|
env | +optional, environmental covariate |
+
generateTestDensityMultiNormal
+ testDensityBanana
+x = c(1,2) +y = testLinearModel(x) +plot(y)
Trace plot for MCMC class
+ + +tracePlot(sampler, thin = "auto", ...)+ +
sampler | +an object of class MCMC sampler |
+
---|---|
thin | +determines the thinning intervall of the chain |
+
... | +additional parameters to pass on to the |
+
marginalPlot
+ plotTimeSeries
+ correlationPlot
+# set up and run the MCMC +ll <- function(x) sum(dnorm(x, log = TRUE)) +setup <- createBayesianSetup(likelihood = ll, lower = c(-10, -10), upper = c(10,10)) +settings <- list(iterations = 2000) +out <- runMCMC(bayesianSetup = setup, settings = settings, sampler = "Metropolis")#> BT runMCMC: trying to find optimal start and covariance values#>#> Running Metropolis-MCMC, chain 1 iteration 100 of 2000 . Current logp: -8.714023 Please wait! Running Metropolis-MCMC, chain 1 iteration 200 of 2000 . Current logp: -8.040798 Please wait! Running Metropolis-MCMC, chain 1 iteration 300 of 2000 . Current logp: -7.858331 Please wait! Running Metropolis-MCMC, chain 1 iteration 400 of 2000 . Current logp: -10.72035 Please wait! Running Metropolis-MCMC, chain 1 iteration 500 of 2000 . Current logp: -10.18255 Please wait! Running Metropolis-MCMC, chain 1 iteration 600 of 2000 . Current logp: -8.271992 Please wait! Running Metropolis-MCMC, chain 1 iteration 700 of 2000 . Current logp: -8.035055 Please wait! Running Metropolis-MCMC, chain 1 iteration 800 of 2000 . Current logp: -9.183765 Please wait! Running Metropolis-MCMC, chain 1 iteration 900 of 2000 . Current logp: -8.860137 Please wait! Running Metropolis-MCMC, chain 1 iteration 1000 of 2000 . Current logp: -8.636479 Please wait! Running Metropolis-MCMC, chain 1 iteration 1100 of 2000 . Current logp: -9.460943 Please wait! Running Metropolis-MCMC, chain 1 iteration 1200 of 2000 . Current logp: -8.746797 Please wait! Running Metropolis-MCMC, chain 1 iteration 1300 of 2000 . Current logp: -9.188869 Please wait! Running Metropolis-MCMC, chain 1 iteration 1400 of 2000 . Current logp: -9.092286 Please wait! Running Metropolis-MCMC, chain 1 iteration 1500 of 2000 . Current logp: -9.186178 Please wait! Running Metropolis-MCMC, chain 1 iteration 1600 of 2000 . Current logp: -9.467755 Please wait! Running Metropolis-MCMC, chain 1 iteration 1700 of 2000 . Current logp: -8.659377 Please wait! Running Metropolis-MCMC, chain 1 iteration 1800 of 2000 . Current logp: -8.451637 Please wait! Running Metropolis-MCMC, chain 1 iteration 1900 of 2000 . Current logp: -14.70914 Please wait! Running Metropolis-MCMC, chain 1 iteration 2000 of 2000 . Current logp: -7.939853 Please wait!#>+# plot the trace +tracePlot(sampler = out, thin = 10)tracePlot(sampler = out, thin = 50)
To update settings of an existing proposal genenerator
+ + +updateProposalGenerator(proposal, chain = NULL, message = F, eps = 1e-10, + manualScaleAdjustment = 1)+ +
proposal | +an object of class proposalGenerator |
+
---|---|
chain | +a chain to create the covariance matrix from (optional) |
+
message | +whether to print an updating message |
+
eps | +numeric tolerance for covariance |
+
manualScaleAdjustment | +optional adjustment for the covariance scale (multiplicative) |
+
The this function can be applied in 2 ways 1) update the covariance given an MCMC chain, and 2) update the proposal generator after parameters have been changed
+ + +Function to plot classic violin plots, as well as "half violin plots" (density plots).
+ + +violinPlot(x, at, .range = 1, add = FALSE, horizontal = TRUE, + which = "both", relToAt = "above", plotQBox = TRUE, plotMed = TRUE, + col = "orangered", border = "black", lwd = par("lwd"), + colQBox = "black", borderQBox = "black", colMed = "white", + pchMed = 19, res = 500, main = "")+ +
x | +vector of values to plot |
+
---|---|
at | +position of the plot when add is TRUE |
+
.range | +maximum height if horizontal or width if vertical of the plot when add is set to FALSE |
+
add | +logical, determining whether the plot should be added to an existing plot window |
+
horizontal | +logical, determining whether the plot should be horizontal, if FALSE the plot will be vertical |
+
which | +string, either "both" for a classic violing plot, or "top" or "bottom" to plot only the half violin |
+
relToAt | +string, one of "centered", "above", "below", "left", "right". Determining the relative position to at. |
+
plotQBox | +logical, determining whether the quantile box should be plotted |
+
plotMed | +logical, determining whether the median should be plotted |
+
col | +color of the violin |
+
border | +color of the border of the violin |
+
lwd | +line width if the border of the violin |
+
colQBox | +color of quantile box |
+
borderQBox | +color of the border of the quantile box |
+
colMed | +color of the median point |
+
pchMed | +pch for median point |
+
res | +"resolution" of the violin. Determining how many descrete points should be used to calculate the density kernel. |
+
main | +header text, only applicable, if 'add' is FALSE |
+