Quick start
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
Installing, loading and citing the package
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.1.
## 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.1},
## url = {https://github.com/florianhartig/BayesianTools},
## }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 Bayesian Setup
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
Running MCMC and SMC functions
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
- A bayesianSetup (alternatively, the log target function)
- The sampler name
- A list with settings - if a parameter is not provided, the default will be used
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)Summarizing outputs
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 codesummary(out)## # # # # # # # # # # # # # # # # # # # # # # # # #
## ## MCMC chain summary ##
## # # # # # # # # # # # # # # # # # # # # # # # # #
##
## # MCMC sampler: Metropolis
## # Nr. Chains: 1
## # Iterations per chain: 10000
## # Rejection rate: 0.686
## # Effective sample size: 932
## # Runtime: 3.378 sec.
##
## # Parameter 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.000and plottted with several plot functions
plot(out) # plot internally calls tracePlot(out)
correlationPlot(out)
marginalPlot(out)## Package 'sm', version 2.2-5.4: type help(sm) for summary information
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).
marginalLikelihood(out)## $marginalLikelihod
## [1] -8.96615
##
## $ln.lik.star
## [1] -2.756817
##
## $ln.pi.star
## [1] -8.987197
##
## $ln.pi.hat
## [1] -2.777863
##
## $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.524434MAP(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.987197You 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 codesummary(out)## # # # # # # # # # # # # # # # # # # # # # # # # #
## ## MCMC chain summary ##
## # # # # # # # # # # # # # # # # # # # # # # # # #
##
## # MCMC sampler: Metropolis
## # Nr. Chains: 3
## # Iterations per chain: 1000
## # Rejection rate: 0.668
## # Effective sample size: 291
## # Runtime: 0.943 sec.
##
## # Parameter MAP 2.5% median 97.5%
## # par 1 : 0.000 -1.829 0.074 1.883
## # par 2 : 0.000 -2.037 0.024 1.935
## # par 3 : 0.000 -2.105 0.063 2.025
##
## ## DIC: 11.579
## ## Convergence
## Gelman Rubin multivariate psrf: 1.003
##
## ## Correlations
## par 1 par 2 par 3
## par 1 1.000 -0.012 0.019
## par 2 -0.012 1.000 0.096
## par 3 0.019 0.096 1.000For 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.01
## par 2 1.01 1.01
## par 3 1.00 1.00
##
## Multivariate psrf
##
## 1Which sampler to choose?
The BT package provides a large class of different MCMC samplers, and it depends on the particular application which is most suitable.
In the absence of further information, we currently recommend the DEzs sampler. This is also the default in the runMCMC function.
