diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 0000000..8dd2c20 --- /dev/null +++ b/.gitattributes @@ -0,0 +1,23 @@ +# Set the default behavior, in case people don't have core.autocrlf set. +* text=auto + +# Explicitly declare text files you want to always be normalized and converted +# to native line endings on checkout. +*.Rd text +*.Rmd text +*.R text +*.r text +*.Rproj text +*.[Rr]md linguist-detectable +*.yml text + + +# Denote all files that are truly binary and should not be modified. +*.Rdata binary +*.RData binary +*.rda binary +*.rdb binary +*.rds binary +*.png binary +*.jpg binary +*.Rdx binary \ No newline at end of file diff --git a/BayesianTools/CRAN.md b/BayesianTools/CRAN.md index b780705..ba91e25 100644 --- a/BayesianTools/CRAN.md +++ b/BayesianTools/CRAN.md @@ -4,15 +4,15 @@ This file contains CRAN submission info / communication ### Submission 1 -This is minor update should re-instate BT to CRAN, after it was removed due to HTML validation problems connected to the recent switch to HTML5 for documentation pages in R 4.2.0, see also https://github.com/florianhartig/BayesianTools/issues/240. +This is a minor update that should bring BT back to CRAN after it was removed due to HTML validation problems associated with the recent switch to HTML5 for documentation pages in R 4.2.0, see also https://github.com/florianhartig/BayesianTools/issues/240. -I had tried to upload a fix to CRAN, which, however, was not accepted due to detritus problems under Windows dev, which, after a long search, turned up to be caused by some stray parallel notes https://github.com/florianhartig/BayesianTools/issues/244 +I had tried to upload a fix to CRAN, but it was not accepted due to detritus problems under Windows dev, which after a long search turned out to be caused by some stray parallel notes https://github.com/florianhartig/BayesianTools/issues/244. -In any case, I hope this will re-instate the BT package to CRAN, sorry for the delay in fixing this. +In any case, I hope this will get the BT package back into CRAN, sorry for the delay in fixing this. See NEWS for other changes. -This release was tested without apparent problems under +This release has been tested without apparent problems on * local MAC OS 13.1, R 4.2.1 Apple * http://win-builder.r-project.org/ - oldrelease / devel / release diff --git a/BayesianTools/DESCRIPTION b/BayesianTools/DESCRIPTION index 0e5dacf..eb4a7ed 100644 --- a/BayesianTools/DESCRIPTION +++ b/BayesianTools/DESCRIPTION @@ -3,7 +3,7 @@ Title: General-Purpose MCMC and SMC Samplers and Tools for Bayesian Statistics Version: 0.1.8 Date: 2023-01-30 Authors@R: c(person("Florian", "Hartig", email = "florian.hartig@biologie.uni-regensburg.de", role = c("aut", "cre"), comment=c(ORCID="0000-0002-6255-9059")), person("Francesco", "Minunno", role = c("aut")), person("Stefan", " Paul", role = c("aut") ), person("David", "Cameron", role = "ctb"), person("Tankred", "Ott", role = "ctb"), person("Maximilian", "Pichler", role = "ctb")) -Description: General-purpose MCMC and SMC samplers, as well as plot and +Description: General-purpose MCMC and SMC samplers, as well as plots and diagnostic functions for Bayesian statistics, with a particular focus on calibrating complex system models. Implemented samplers include various Metropolis MCMC variants (including adaptive and/or delayed rejection MH), the diff --git a/BayesianTools/R/BayesianSetupGenerateParallel.R b/BayesianTools/R/BayesianSetupGenerateParallel.R index 4a5cb5a..766538f 100644 --- a/BayesianTools/R/BayesianSetupGenerateParallel.R +++ b/BayesianTools/R/BayesianSetupGenerateParallel.R @@ -1,18 +1,13 @@ #' Factory to generate a parallel executor of an existing function -#' #' @author Florian Hartig #' @param fun function to be changed to parallel execution -#' @param 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 -#' @param parallelOptions list containing three lists. First "packages" determines the R packages necessary to run the likelihood function. Second "variables" the objects in the global environment needed to run the likelihood function and third "dlls" the DLLs needed to run the likelihood function (see Details). -#' @note Can also be used to make functions compatible with library sensitivity -#' @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()). +#' @param parallel should a parallel R cluster be used? If set to T, the operating system will automatically detect the available cores and n-1 of the available n cores will be used. Alternatively, you can manually set the number of cores to be used +#' @param parallelOptions a list containing three lists. \itemize{\item First, "packages": determines the R packages required to run the likelihood function. \item Second, "variables": the objects in the global environment needed to run the likelihood function. \item Third, "dlls": the DLLs needed to run the likelihood function (see Details).} +#' @note can be used to make functions compatible with library sensitivity +#' @details For parallelization, if option T is selected, an automatic parallelization is tried via R. Alternatively, "external" can be selected on the assumption that the likelihood has already been parallelized. In the latter case, a matrix with parameters as columns must be accepted. You can also specify which packages, objects and DLLs are exported to the cluster. By default, a copy of your workspace is exported, but depending on your workspace, this can be inefficient. As an alternative, you can specify the environments and packages in the likelihood function (e.g. BayesianTools::VSEM() instead of VSEM()). #' @export #' @example /inst/examples/generateParallelExecuter.R - +#' generateParallelExecuter <- function(fun, parallel = F, parallelOptions = list(variables = "all", packages = "all", dlls = NULL)){ if (parallel == F){ @@ -27,7 +22,6 @@ generateParallelExecuter <- function(fun, parallel = F, parallelOptions = list(v #library(foreach) #library(iterators) # library(parallel) - if (parallel == T | parallel == "auto"){ cores <- parallel::detectCores() - 1 } else if (is.numeric(parallel)){ diff --git a/BayesianTools/R/BayesianTools.R b/BayesianTools/R/BayesianTools.R index 3b79b19..701f39e 100644 --- a/BayesianTools/R/BayesianTools.R +++ b/BayesianTools/R/BayesianTools.R @@ -1,5 +1,6 @@ #' @title BayesianTools #' @name BayesianTools +#' @aliases BayesianTools-package #' @docType package #' @useDynLib BayesianTools, .registration = TRUE #' @description A package with general-purpose MCMC and SMC samplers, as well as plots and diagnostic functions for Bayesian statistics diff --git a/BayesianTools/R/classBayesianSetup.R b/BayesianTools/R/classBayesianSetup.R index 5441e59..0fae73b 100644 --- a/BayesianTools/R/classBayesianSetup.R +++ b/BayesianTools/R/classBayesianSetup.R @@ -1,3 +1,4 @@ + #' Creates a standardized collection of prior, likelihood and posterior functions, including error checks etc. #' @author Florian Hartig, Tankred Ott #' @param likelihood log likelihood density function diff --git a/BayesianTools/R/classMcmcSampler.R b/BayesianTools/R/classMcmcSampler.R index 1f34d7e..8d9fc3e 100644 --- a/BayesianTools/R/classMcmcSampler.R +++ b/BayesianTools/R/classMcmcSampler.R @@ -116,7 +116,7 @@ getSample.mcmcSampler <- function(sampler, parametersOnly = T, coda = F, start = #' @method summary mcmcSampler #' @author Stefan Paul #' @export -summary.mcmcSampler <- function(object, ...){ +summary.mcmcSampler <- function(object, printCorrelation = "auto", ...){ #codaChain = getSample(sampler, parametersOnly = parametersOnly, coda = T, ...) #summary(codaChain) #rejectionRate(sampler$codaChain) @@ -133,8 +133,7 @@ summary.mcmcSampler <- function(object, ...){ mcmcsampler <- sampler$settings$sampler runtime <- sampler$settings$runtime[3] - correlations <- round(cor(getSample(sampler)),3) - + chain <- getSample(sampler, parametersOnly = T, coda = T, ...) # chain <- getSample(sampler, parametersOnly = T, coda = T) if("mcmc.list" %in% class(chain)){ @@ -203,8 +202,11 @@ summary.mcmcSampler <- function(object, ...){ try(cat("## DIC: ", round(DInf$DIC,3), "\n"), silent = TRUE) cat("## Convergence" ,"\n", "Gelman Rubin multivariate psrf: ", conv, "\n","\n") - cat("## Correlations", "\n") - print(correlations) + if(printCorrelation == TRUE){ + correlations <- round(cor(getSample(sampler)),3) + cat("## Correlations", "\n") + print(correlations) + } } #' @author Florian Hartig diff --git a/BayesianTools/R/mcmcRun.R b/BayesianTools/R/mcmcRun.R index 91f7b9b..5258cef 100644 --- a/BayesianTools/R/mcmcRun.R +++ b/BayesianTools/R/mcmcRun.R @@ -263,7 +263,30 @@ runMCMC <- function(bayesianSetup , sampler = "DEzs", settings = NULL){ #' @param settings optional list with parameters that will be used instead of the defaults #' @param sampler one of the samplers in \code{\link{runMCMC}} #' @param check logical determines whether parameters should be checked for consistency -#' @details see \code{\link{runMCMC}} +#' @details +#' +#' The following settings can be used for all MCMCs: +#' +#' startValue (no default) start values for the MCMC. Note that DE family samplers require a matrix of #' start values. If startvalues are not provided, they are sampled from the prior. +#' +#' iterations (10000) the MCMC iterations +#' +#' burnin (0) burnin +#' +#' thin (1) thinning while sampling +#' +#' consoleUpdates (100) update frequency for console updates +#' +#' parallel (NULL) whether parallelization is to be used +#' +#' message (TRUE) if progress messages are to be printed +#' +#' nrChains (1) the number of independent MCMC chains to be run. Note that this is not controlling the #' internal number of chains in population MCMCs such as DE, so if you run nrChains = 3 with a DEzs #' #' startValue that is a 4xparameter matrix (= 4 internal chains), you will run independent DEzs runs #' #' with 4 internal chains each. +#' +#' For more details, see \code{\link{runMCMC}} +#' +#' +#' @example inst/examples/mcmcRun.R #' @export applySettingsDefault<-function(settings=NULL, sampler = "DEzs", check = FALSE){ diff --git a/BayesianTools/R/plotMarginals.R b/BayesianTools/R/plotMarginals.R index b962779..ca829a4 100644 --- a/BayesianTools/R/plotMarginals.R +++ b/BayesianTools/R/plotMarginals.R @@ -170,7 +170,7 @@ marginalPlotDensity <- function(posteriorMat, priorMat = NULL, xrange = NULL, co mtext('Marginal parameter uncertainty', outer = TRUE, cex = 1.5) } else { - mfrow <- if (nPar < 16) getPanels(nPar) else c(4,4) + mfrow <- getPanels(nPar) op <- par(mfrow = mfrow, mar = c(4.5, 4, 5, 3), oma=c(3, 1.5, 2, 0), xpd=TRUE) on.exit(par(op)) @@ -282,7 +282,7 @@ marginalPlotViolin <- function(posteriorMat, priorMat = NULL, xrange = NULL, col mtext('Marginal parameter uncertainty', outer = TRUE, cex = 1.5) } else { - mfrow <- if (nPar < 16) getPanels(nPar) else c(4,4) + mfrow <- getPanels(nPar) op <- par(mfrow = mfrow, mar = c(4.5, 4.5, 5, 3), oma=c(3, 0, 2, 0), xpd=TRUE) diff --git a/BayesianTools/inst/examples/mcmcRun.R b/BayesianTools/inst/examples/mcmcRun.R index 7e5ea09..029e670 100644 --- a/BayesianTools/inst/examples/mcmcRun.R +++ b/BayesianTools/inst/examples/mcmcRun.R @@ -5,8 +5,9 @@ ll <- generateTestDensityMultiNormal(sigma = "no correlation") ## 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) +## Finally we can run the sampler. To get possible settings +## for a sampler, see help or run applySettingsDefault(sampler = "Metropolis") +settings = list(iterations = 1000, adapt = FALSE) # out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings) ## out is of class bayesianOutput. There are various standard functions @@ -21,3 +22,5 @@ summary(out) # for plotting and analysis codaObject = getSample(out, start = 500, coda = TRUE) + + diff --git a/BayesianTools/man/BayesianTools.Rd b/BayesianTools/man/BayesianTools.Rd index 2c5f317..9ab049a 100644 --- a/BayesianTools/man/BayesianTools.Rd +++ b/BayesianTools/man/BayesianTools.Rd @@ -3,6 +3,7 @@ \docType{package} \name{BayesianTools} \alias{BayesianTools} +\alias{BayesianTools-package} \title{BayesianTools} \description{ A package with general-purpose MCMC and SMC samplers, as well as plots and diagnostic functions for Bayesian statistics diff --git a/BayesianTools/man/applySettingsDefault.Rd b/BayesianTools/man/applySettingsDefault.Rd index 5e7290e..f1abf64 100644 --- a/BayesianTools/man/applySettingsDefault.Rd +++ b/BayesianTools/man/applySettingsDefault.Rd @@ -17,7 +17,53 @@ applySettingsDefault(settings = NULL, sampler = "DEzs", check = FALSE) Provides the default settings for the different samplers in runMCMC } \details{ -see \code{\link{runMCMC}} +The following settings can be used for all MCMCs: + +startValue (no default) start values for the MCMC. Note that DE family samplers require a matrix of #' start values. If startvalues are not provided, they are sampled from the prior. + +iterations (10000) the MCMC iterations + +burnin (0) burnin + +thin (1) thinning while sampling + +consoleUpdates (100) update frequency for console updates + +parallel (NULL) whether parallelization is to be used + +message (TRUE) if progress messages are to be printed + +nrChains (1) the number of independent MCMC chains to be run. Note that this is not controlling the #' internal number of chains in population MCMCs such as DE, so if you run nrChains = 3 with a DEzs #' #' startValue that is a 4xparameter matrix (= 4 internal chains), you will run independent DEzs runs #' #' with 4 internal chains each. + +For more details, see \code{\link{runMCMC}} +} +\examples{ +## 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. To get possible settings +## for a sampler, see help or run applySettingsDefault(sampler = "Metropolis") +settings = list(iterations = 1000, adapt = FALSE) # +out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings) + +## out is of class bayesianOutput. There are various standard functions +# implemented for this output + +plot(out) +correlationPlot(out) +marginalPlot(out) +summary(out) + +## 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) + + } \author{ Florian Hartig diff --git a/BayesianTools/man/generateParallelExecuter.Rd b/BayesianTools/man/generateParallelExecuter.Rd index 82b7a37..c41ae3f 100644 --- a/BayesianTools/man/generateParallelExecuter.Rd +++ b/BayesianTools/man/generateParallelExecuter.Rd @@ -13,22 +13,18 @@ generateParallelExecuter( \arguments{ \item{fun}{function to be changed to parallel execution} -\item{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} +\item{parallel}{should a parallel R cluster be used? If set to T, the operating system will automatically detect the available cores and n-1 of the available n cores will be used. Alternatively, you can manually set the number of cores to be used} -\item{parallelOptions}{list containing three lists. First "packages" determines the R packages necessary to run the likelihood function. Second "variables" the objects in the global environment needed to run the likelihood function and third "dlls" the DLLs needed to run the likelihood function (see Details).} +\item{parallelOptions}{a list containing three lists. \itemize{\item First, "packages": determines the R packages required to run the likelihood function. \item Second, "variables": the objects in the global environment needed to run the likelihood function. \item Third, "dlls": the DLLs needed to run the likelihood function (see Details).}} } \description{ Factory to generate a parallel executor of an existing function } \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()). +For parallelization, if option T is selected, an automatic parallelization is tried via R. Alternatively, "external" can be selected on the assumption that the likelihood has already been parallelized. In the latter case, a matrix with parameters as columns must be accepted. You can also specify which packages, objects and DLLs are exported to the cluster. By default, a copy of your workspace is exported, but depending on your workspace, this can be inefficient. As an alternative, you can specify the environments and packages in the likelihood function (e.g. BayesianTools::VSEM() instead of VSEM()). } \note{ -Can also be used to make functions compatible with library sensitivity +can be used to make functions compatible with library sensitivity } \examples{ diff --git a/BayesianTools/man/runMCMC.Rd b/BayesianTools/man/runMCMC.Rd index 20b6687..eacee0f 100644 --- a/BayesianTools/man/runMCMC.Rd +++ b/BayesianTools/man/runMCMC.Rd @@ -56,8 +56,9 @@ ll <- generateTestDensityMultiNormal(sigma = "no correlation") ## 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) +## Finally we can run the sampler. To get possible settings +## for a sampler, see help or run applySettingsDefault(sampler = "Metropolis") +settings = list(iterations = 1000, adapt = FALSE) # out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings) ## out is of class bayesianOutput. There are various standard functions @@ -72,6 +73,8 @@ summary(out) # for plotting and analysis codaObject = getSample(out, start = 500, coda = TRUE) + + } \seealso{ \code{\link{createBayesianSetup}} diff --git a/BayesianTools/vignettes/BayesianTools.Rmd b/BayesianTools/vignettes/BayesianTools.Rmd index 42484d1..8ad9f5d 100644 --- a/BayesianTools/vignettes/BayesianTools.Rmd +++ b/BayesianTools/vignettes/BayesianTools.Rmd @@ -5,9 +5,13 @@ output: toc: true vignette: > %\VignetteIndexEntry{Manual for the BayesianTools R package} - %\VignetteEngine{knitr::rmarkdown} \usepackage[utf8]{inputenc} + %\VignetteEngine{knitr::rmarkdown} 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." +editor_options: + markdown: + wrap: 72 + chunk_output_type: console --- ```{r global_options, include=FALSE} @@ -20,9 +24,9 @@ set.seed(123) # 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 +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 following sections. -## Installing, loading and citing the package +### Install, load and cite the package If you haven't installed the package yet, either run @@ -39,75 +43,109 @@ library(BayesianTools) citation("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! +Note: BayesianTools calls a number of secondary packages. Of particular importance is `coda`, which is used in a number of plots and summarystatistics. If you use a lot of summary statistics and diagnostic 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) +Pro Tip: If you are running a stochastic algorithm such as MCMC, you should always set or record your random seed to make your results reproducible (otherwise the results will change slightly each time you run the code). ```{r} set.seed(123) ``` -In a real application, to ensure reproducibility, it would also be useful to record the session info, +In a real application, recording the session info would be helpful in ensuring reproducibility. ```{r, eval = F} sessionInfo() ``` -which lists the version number of R and all loaded packages. +This session information includes 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. +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. + +The `createBayesianSetup` function generates a `BayesianSetup` object. The function requires a log-likelihood and a log-prior (optional) as arguments. It then automatically creates the posterior and variousn convenience functions for the samplers. -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 -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). + 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. +If no prior information is provided, an unbounded flat prior is generated. If no explicit prior is specified, but lower and upper values are given, a standard uniform prior with the respective bounds is created, including the option to sample from this prior, which is useful for SMC and obtaining initial values. This option is used in the following example, which creates a multivariate normal likelihood density and a uniform prior for 3 parameters. ```{r} 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. +A more detailed description of the `BayesianSetup` will follow below. -**Hint:** for an example how to run this steps for dynamic ecological model, see ?VSEM +**Hint:** For an example of how to run these steps for a dynamic ecological model, see `?VSEM`. -## Running MCMC and SMC functions +## Run 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 +After setup, you may need to run a calibration. The `runMCMC` function serves as the main wrapper for all other implemented `MCMC`/`SMC` functions. It always requires 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 +- `bayesianSetup` (alternatively, the log target function) +- `sampler` name +- list with `settings` for each sampler - if `settings` is not + specified, the default value will be applied -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 later reference on MCMC samplers. This is how we would call this sampler with default settings +For example, choosing the `sampler` name "Metropolis" calls a versatile Metropolis-type MCMC with options for covariance adjustment, delayed rejection, tempering and Metropolis-within-Gibbs sampling. For details, see the later reference on MCMC samplers. When in doubt about the choice of MCMC sampler, we recommend using the default "DEzs". ```{r} -iter = 10000 +iter = 1000 settings = list(iterations = iter, message = FALSE) out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings) ``` -#### Summarizing outputs +## Convergence checks for MCMCs -All samplers can be plotted and summarized via the console with the standard print, and summary commands +Before interpreting the results, MCMCs should be checked for convergence. We recommend the Gelman-Rubin convergence diagnostics,which are the standard used in most publications. The Gelman-Rubin diagnostics requires running multiple MCMCs (we recommend 3-5). + +For all samplers, you can conveniently perform multiple runs using the `nrChains` argument. Alternatively, for runtime reasons, you can combine the results of three independent `runMCMC` runs with `nrChains = 1` using `combineChains` + +```{r, echo = T} +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 you to do everything you can do with an `mcmcSampler` object (sometimes with slightly different output). + +Basic convergence is checked via so-called trace plots, which show the 3 MCMC chains in different colors. ```{r} -print(out) -summary(out) +plot(out) ``` -and plotted 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"). +The trace plot is examined for major problems (chains look different, systematic trends) and to decide on the burn-in, i.e. the point at which the MCMC has reached the sampling range (i.e. the chains no longer systematically go in one direction). Here, they have basically reached this point immediately, so we could set the burn-in to 0, but I choose 100, i.e. discard the first 100 samples in all further plots. + +If the trace plots look good, we can look at the Gelman-Rubin convergence diagnostics. Note that you have to discard the burn-in. + +```{r} +gelmanDiagnostics(out, plot = F, start = 100) +``` + +Usually, a value \< 1.05 for each parameter and a msrf \< 1.1 of 1.2 are considered sufficient for convergence. + +### Summarize the outputs + +If we are happy with the convergence, we can plot and summarize all `sampler`s from the console with the standard `print` and `summary` commands. + +```{r} +print(out, start = 100) +summary(out, start = 100) +``` + +You can also use built-in plot functions from the package for visualization. The `marginalPlot` can be either plotted as histograms with density overlay (default setting) or as a violin plot (see "?marginalPlot"). ```{r} -plot(out) # plot internally calls tracePlot(out) correlationPlot(out) marginalPlot(out, prior = TRUE) ``` -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"). +Additional functions that may be used on all `sampler`s are model selection scores, including the Deviance Information Criterion (DIC) and the marginal likelihood (see later section for details on calculating the Bayes factor), alongside the Maximum A Posteriori (MAP) value. A set of methods are available for calculation of marginal likelihood (refer to "?marginalLikelihood"). ```{r} marginalLikelihood(out) @@ -115,13 +153,13 @@ DIC(out) MAP(out) ``` -You can extract (a part of) the sampled parameter values by +To extract part of the sampled parameter values, you can use the following process: ```{r, eval = F} getSample(out, start = 100, end = NULL, thin = 5, whichParameters = 1:2) ``` -For all samplers, you can conveniently perform multiple runs via the nrChains argument +### Which sampler to choose? ```{r, echo = T} iter = 1000 @@ -152,9 +190,9 @@ gelmanDiagnostics(out, plot = T) #### Which 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. +The BT package provides a large class of different MCMC samplers, and it depends on the particular application which one is most suitable. -In the absence of further information, we currently recommend the DEzs sampler. This is also the default in the runMCMC function. +In the absence of further information, we currently recommend the `DEz`s sampler. This is also the default in the `runMCMC` function. # BayesianSetup Reference @@ -166,38 +204,38 @@ 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. +See parallelization options below. We will use a simple 3-d multivariate normal density for the demonstration. ```{r} ll <- generateTestDensityMultiNormal(sigma = "no correlation") bayesianSetup = createBayesianSetup(likelihood = ll, lower = rep(-10, 3), upper = rep(10, 3)) ``` -#### Parallelization of the likelihood evaluations +### Parallelization of the likelihood evaluations -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 +Likelihoods are often costly to compute. If this is the case for you, you should think about parallelization possibilities. The 'createBayesianSetup' function has an input variable 'parallel', with the following options - F / FALSE means no parallelization should be used -- T / TRUE means that automatic parallelization options from R are used (careful: this will not work if your likelihood writes to file, or uses global variables or functions - see general R help on parallelization) -- "external", assumed that the likelihood is already parallelized. In this case, the function needs to accept a matrix with parameters as columns, and rows as the different model runs you want to evaluate. This is the most likely option to use if you have a complicated setup (file I/O, HPC cluster) that cannot be treated with the standard R parallelization. +- T / TRUE means to use R's automatic parallelization options (note: this will not work if your likelihood is writing to file, or using global variables or functions - see R's general help on parallelization) +- "external", assumes that the likelihood is already parallelized. In this case, the function needs to take a matrix with parameters as columns, and rows as the different model runs you want to evaluate. This is the most likely option to use if you have a complicated setup (file I/O, HPC cluster) that cannot be handled with the standard R parallelization. -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. +Algorithms in the `BayesianTools` package can take advantage of parallelization if this option is specified in the `BayesianSetup`. Note that parallelization is currently 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 +Here are some more details about the parallelization. #### 1. In-build 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. +In-built parallelizing is the easiest way to use parallel computing. The "parallel" argument allows you to select the number of cores to use for parallelization. Alternatively, TRUE or "auto" will use all availablecores except one. Now the proposals are evaluated in parallel. Technically, the built-in parallelization uses an R cluster to evaluate the posterior density function. The input to the parallel function is a matrix where each column represents a parameter and each row represents a proposal. This allows the proposals to be evaluated in parallel. For `sampler`s that evaluate only one proposal at a time (namely the Metropolis-based algorithms and DE/DREAM without the `zs` extension), parallelization cannot be used. #### 2. External parallelization -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: +The second option is to use external parallelization. Here, parallelization is attempted in the user-defined likelihood function. To use external parallelization, the likelihood function must take a matrix of proposals and return a vector of likelihood values. In the proposal matrix, each row represents a proposal, and each column represents a parameter. In addition, you will need to specify the "external" parallelization in the "parallel" argument. In simple terms, using external parallelization involves the following steps ```{r, eval = FALSE} ## Definition of likelihood function likelihood <- function(matrix){ # Calculate likelihood in parallel - # Return vector of likelihood valus + # Return vector of likelihood values } ## Create Bayesian Setup @@ -209,13 +247,13 @@ runMCMC(BS, sampler = "SMC", ...) #### 3. Multi-core and cluster calculations -If you want to run your calculations on a cluster there are several ways to achieve it. +If you want to run your calculations on a cluster there are several ways to do it. -In the first case you want to parallize n internal (not overall chains) on n cores. The argument "parallel = T" in "createBayesianSetup" allows only at most parallelization on 3 cores for the SMC, DEzs and DreamsSamplers. But by setting "parallel = n" to n cores in the "createBayesianSetup", the internal chains of DEzs and DREAMzs will be parallelized on n cores. This works only for the DEzs and DREAMzs samplers. +In the first case, you want to parallelize n internal (not total chains) on n cores. The argument "parallel = T" in "createBayesianSetup" only allows parallelization on a maximum of 3 cores for the `SMC`, `DEzs` and `DreamsSamplers`. But by setting "parallel = n" in "createBayesianSetup" to `n` cores, the internal chains of `DEzs` and `DREAMzs` will be parallelized on `n` cores. This only works for `DEzs` and `DREAMzs` samplers. ```{r, eval = FALSE} ## n = Number of cores -n=2 +n = 2 x <- c(1:10) likelihood <- function(param) return(sum(dnorm(x, mean = param, log = T))) bayesianSetup <- createBayesianSetup(likelihood, parallel = n, lower = -5, upper = 5) @@ -224,7 +262,7 @@ bayesianSetup <- createBayesianSetup(likelihood, parallel = n, lower = -5, upper out <- runMCMC(bayesianSetup, settings = list(iterations = 1000)) ``` -In the second case you want to parallize n internal chains on n cores with a external parallilzed likelihood function. Unlike the previous case, that way DEzs, DREAMzs, and SMC samplers can be parallelized. +In the second case, you want to parallelize n internal chains on `n` cores with an external parallelized likelihood function. Unlike the previous case, `DEzs`, `DREAMzs`, and `SMC` samplers can be parallelized this way. ```{r, eval = FALSE} ### Create cluster with n cores @@ -250,7 +288,8 @@ settings = list(iterations = 100, nrChains = 1, startValue = bayesianSetup$prior out <- runMCMC(bayesianSetup, settings, sampler = "DEzs") ``` -In a another case your likelihood requires a parallized model. Start your cluster and export your model, the required libraries, and dlls. Now you can start your calculations with the argument "parallel = external" in createBayesianSetup. +In another case, your likelihood requires a parallelized model. Start your cluster and export your model, the required libraries, and `dll`s. Now you can start your calculations with the argument "parallel = +external" in `createBayesianSetup`. ```{r, eval = FALSE} ### Create cluster with n cores @@ -274,7 +313,7 @@ out <- runMCMC(bayesianSetup, settings) ``` -In the last case you can parallize over whole chain calculations. However, here the likelihood itself will not be parallelized. Each chain will be run on one core and the likelihood will be calculated on that core. +In the last case, you can parallelize over whole chain calculations. However, the likelihood itself is not parallelized. Each chain is run on one core, and the likelihood is calculated on that core. ```{r, eval = FALSE} ### Definition of likelihood function @@ -296,18 +335,19 @@ out <- parallel::parLapply(cl, 1:n, fun = function(X, bayesianSetup, settings) r out <- createMcmcSamplerList(out) ``` -\*\* 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. \*\* +\*\* Note: 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 information from the parallel cores). For models with low computational cost, this procedure may take more time than the actual evaluation of the likelihood. If in doubt, do a small run-time comparison before starting your large sampling. \*\* ## Reference on creating priors -The prior in the BayesianSetup consists of four parts +The prior in the `BayesianSetup` consists of four parts - A log density function -- An (optional) sampling function (must be a function without parameters, that returns a draw from the prior) +- An (optional) sampling function (must be a function without + parameters, that returns a draw from the prior) - lower / upper boundaries - Additional info - best values, names of the parameters, ... -These information can be passed by first creating an extra object, via createPrior, or through the the createBayesianSetup function. +This information can be passed by first creating an extra object, via `createPrior`, or via the `createBayesianSetup` function. #### Creating priors @@ -315,13 +355,13 @@ You have 5 options to create a prior - Do not set a prior - in this case, an infinite prior will be created - Set min/max values - a bounded flat prior and the corresponding sampling function will be created -- Use one of the pre-definded priors, see ?createPrior for a list. One of the options here is to use a previous MCMC output as new prior. Pre-defined priors will usually come with a sampling function -- Use a user-defined prior, see ?createPrior +- Use one of the predefinded priors, see `?createPrior` for a list. One of the options here is to use a previous MCMC output as the new prior. Predefined priors will usually come with a sampling function +- Use a user-defined prior, see `?createPrior` - Create a prior from a previous MCMC sample #### Creating user-defined priors -If creating a user-defined prior, the following information can/should be provided to createPrior: +When creating a user-defined prior, the following information can/should be passed to `createPrior` - A log density function, as a function of a parameter vector x, same syntax as the likelihood - Additionally, you should consider providing a function that samples from the prior, because many samplers (SMC, DE, DREAM) can make use of this function for initial conditions. If you use one of the pre-defined priors, the sampling function is already implemented @@ -330,7 +370,7 @@ If creating a user-defined prior, the following information can/should be provid #### Creating a prior from a previous MCMC sample -The following example from the help shows how this works +The following example from the help file illustrates the process ```{r} # Create a BayesianSetup @@ -352,11 +392,13 @@ out <- runMCMC(bayesianSetup = bayesianSetup, settings = settings) ## The runMCMC() function -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. +The `runMCMC` function is the central function for starting MCMC algorithms in the `BayesianTools` package. It takes a `bayesianSetup`, a choice of `sampler` (default is `DEzs`), and optionally changes to the default settings of the chosen `sampler`. +```{r, message = F} 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. +You may use an optional argument `nrChains`, which is set to the default value of 1 but can be modified if needed. Increasing its value will lead `runMCMC` to execute multiple runs. ```{r, message = F} @@ -392,7 +434,7 @@ newPriorFromPosterior <- createPriorDensity(out2) ## The different MCMC samplers -For convenience we define a number of iterations +For simplicity, we will define a fixed number of iterations. ```{r} iter = 10000 @@ -400,19 +442,20 @@ iter = 10000 ### The Metropolis MCMC class -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 `BayesianTools` package is able to run a large number of Metropolis-Hastings (MH) based algorithms. All of these `sampler`s 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. +The subsequent code provides an overview of the default settings of the MH `sampler`. ```{r} applySettingsDefault(sampler = "Metropolis") ``` -The following examples show how the different settings can be used. As you will see different options can be activated singly or in combination. +Here are some examples of how to apply different settings. Activate individual options or combinations as demonstrated. -#### Standard MH MCMC +### Standard MH MCMC -The following settings will run the standard Metropolis Hastings MCMC. +The following settings 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. @@ -426,7 +469,7 @@ plot(out) #### Standard MH MCMC, prior optimization -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. +Prior to the sampling process, this `sampler` employs an optimization step. The purpose of optimization is to improve the initial values and the covariance of the proposal distribution. ```{r, results = 'hide', eval = F} settings <- list(iterations = iter, adapt = F, DRlevels = 1, gibbsProbabilities = NULL, temperingFunction = NULL, optimize = T, message = FALSE) @@ -436,7 +479,7 @@ plot(out) #### Adaptive MCMC, prior optimization -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. +The adaptive Metropolis sampler (AM) uses the information already obtained during the sampling process to improve (or adapt) the proposal function. The `BayesianTools` package adjusts the covariance of the proposal distribution by utilizing the chain's history. References: Haario, H., E. Saksman, and J. Tamminen (2001). An adaptive metropolis algorithm. Bernoulli , 223-242. @@ -448,7 +491,7 @@ plot(out) #### Standard MCMC, prior optimization, delayed rejection -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. \*\* +Although rejection is an essential step in an MCMC algorithm, it can also mean that the proposal distribution is (locally) poorly 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 greater flexibility of the sampler. In the `BayesianTools` package, the number of delayed rejection steps and the scaling of the proposals can be specified. \*\* Note that the current version supports only two delayed rejection steps. \*\* References: Green, Peter J., and Antonietta Mira. "Delayed rejection in reversible jump Metropolis-Hastings." Biometrika (2001): 1035-1053. @@ -460,7 +503,7 @@ plot(out) #### Adaptive MCMC, prior optimization, delayed rejection -The delayed rejection adaptive Metropolis (DRAM) sampler is merely a combination of the two previous sampler (DR and AM). +The Delayed Rejection Adaptive Metropolis (DRAM) `sampler` is simply a combination of the two previous `sampler`s (DR and AM). References: Haario, Heikki, et al. "DRAM: efficient adaptive MCMC." Statistics and Computing 16.4 (2006): 339-354. @@ -472,9 +515,10 @@ plot(out) #### Standard MCMC, prior optimization, Gibbs updating -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. +To reduce the dimensions of the target function, a Metropolis-within-Gibbs `sampler` can be run with the can be run with the `BayesianTools` package. This means that only a subset of the parameter vector is updated in each iteration. In the example below, at most two (of the three) parameters are updated at each step, and it is twice as likely to vary one as to vary two. -\*\* Note that currently adaptive cannot be mixed with Gibbs updating! \*\* +\*\* Note that currently adaptive cannot be mixed with Gibbs updating! +\*\* ```{r, results = 'hide', eval = F} settings <- list(iterations = iter, adapt = T, DRlevels = 1, gibbsProbabilities = c(1,0.5,0), temperingFunction = NULL, optimize = T, message = FALSE) @@ -484,7 +528,7 @@ plot(out) #### Standard MCMC, prior optimization, gibbs updating, tempering -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. +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 lead to a faster initial scanning of the target function. To incorporate this, a tempering function must be provided by the user. The function describes how to influence the acceptance rate during burn-in. In the example below, an exponential decline approaching 1 (= no influence 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. @@ -499,9 +543,9 @@ plot(out) ### Differential Evolution MCMC -The BT package implements two versions of the differential evolution MCMC. In doubt, you should use the DEzs option. +The BT package implements two differential evolution MCMCs. When in doubt, 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 difference to the Metropolis based algorithms is the creation of the proposal. 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. +The first is the normal DE MCMC, according 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. This sampler runs multiple chains in parallel (but not in the sense of parallel computing). The main difference to the Metropolis based algorithms is the generation of the proposal. In general, all `sampler`s use the current position of the chain and add a step in the parameter space to generate a new proposal. While in the Metropolis based `sampler` this step is usually drawn from a multivariate normal distribution (but any distribution is possible), the `DE` `sampler` uses the current position of two other chains to generate the step for each chain. For successful sampling, at least `2*d` chains, where `d` is the number of parameters, must be run in parallel. ```{r, results = 'hide', eval = F} settings <- list(iterations = iter, message = FALSE) @@ -509,7 +553,7 @@ out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "DE", settings = setting 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. +The second is the Differential Evolution MCMC with Snooker update and sampling from past states, according 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 from the normal `DE` MCMC. First, it uses a snooker update based on a user-defined probability. Second, past states of other chains are taken into account when generating the proposal. These extensions allow fewer chains (i.e., 3 chains are usually sufficient for up to 200 parameters) and parallel computing, since the current position of each chain depends only on the past states of the other chains. ```{r, results = 'hide', eval = F} settings <- list(iterations = iter, message = FALSE) @@ -517,11 +561,14 @@ out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "DEzs", settings = setti plot(out) ``` -### DREAM sampler +## DREAM sampler + +There are two versions of the DREAM `sampler`. First, 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. -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 based on the DE sampler with some significant changes: -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. + 1) More than two chains can be used to generate a proposal. + 2) Randomized subspace sampling can be used to improve efficiency for high-dimensional posteriors. Each dimension is updated with a crossover probability CR. To speed up the exploration of the posterior, DREAM adjusts the distribution of CR values during burn-in to favor large jumps over small ones. 3) Outlier chains can be removed during burn-in. ```{r, results = 'hide', eval = F} settings <- list(iterations = iter, message = FALSE) @@ -529,9 +576,9 @@ out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "DREAM", settings = sett 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. +The second implementation uses the same extension as the DEzs sampler. Namely sampling from past states and a snooker update. Again, this extension allows the use of fewer chains and parallel computing. -Again, in doubt you should prefer "DREAMzs". +Again, if in doubt, you should prefer "DREAMzs". ```{r, results = 'hide', eval = FALSE} settings <- list(iterations = iter, message = FALSE) @@ -539,9 +586,9 @@ out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "DREAMzs", settings = se plot(out) ``` -### T-walk +## T-walk -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. +The t-walk is an 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. The sampler uses two independent points to explore the posterior space. Based on probabilities, four different moves are used to generate proposals for the two points. As with the DE sampler, this procedure does not require tuning of the proposal distribution for efficient sampling in complex posterior distributions. ```{r, eval = F} settings = list(iterations = iter, message = FALSE) @@ -549,14 +596,13 @@ settings = list(iterations = iter, message = FALSE) out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Twalk", settings = settings) ``` -### Convergence checks for MCMCs +## Non-MCMC sampling algorithms -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. +MCMCs sample the posterior space by creating a chain in parameter space. While this allows for "learning" from past steps, it does not allow for running a large number of posteriors in parallel. -```{r, eval = T} -settings <- list(iterations = iter, nrChains = 3, message = FALSE) -out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings) -plot(out) +An alternative to MCMCs are particle filters, also known as 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 + +### Rejection sampling #chain = getSample(out, coda = T) gelmanDiagnostics(out, plot = F) @@ -570,7 +616,7 @@ An alternative to MCMCs are particle filters, aka Sequential Monte-Carlo (SMC) a ### Rejection sampling -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. +The simplest option is to sample a large number of parameters and accept them according to their posterior value. This option can be emulated with the implemented SMC by setting iterations to 1. ```{r, results = 'hide', eval = F} settings <- list(initialParticles = iter, iterations= 1) @@ -580,7 +626,7 @@ plot(out) ### Sequential Monte Carlo (SMC) -The more sophisticated option is using the implemented SMC, which is basically a particle filter that applies several filter steps. +The more sophisticated option is to use the implemented SMC, which is basically a particle filter that applies multiple filter steps. ```{r, results = 'hide', eval = F} settings <- list(initialParticles = iter, iterations= 10) @@ -588,27 +634,27 @@ out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "SMC", settings = settin plot(out) ``` -Note that the use of a number for initialParticles requires that the bayesianSetup includes the possibility to sample from the prior. +Note that using a number for `initialParticles` requires that the `bayesianSetup` includes the option to sample from the prior. # Bayesian model comparison and averaging -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. +There are a number of Bayesian methods for model selection and model comparison. The BT implements three of the most common: DIC, 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. +- For an overview of DIC and WAIC, see 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 Bayes factor relies on the estimation of marginal likelihoods, which is numerically challenging. 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 recommended method is the "Chib" method (Chib and Jeliazkov, 2001), which is based on MCMC samples but performs additional calculation. Although this is the current recommendation, note that there are some numerical issues with this algorithm that may limit its 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 third method is simply sampling from the prior. While in principle unbiased, it converges only for a large number of samples and is therefore numerically inefficient. ## Example -Data linear Regression with quadratic and linear effect +Data linear regression with quadratic and linear effect ```{r} sampleSize = 30 @@ -641,7 +687,7 @@ setUp1 <- createBayesianSetup(likelihood1, lower = c(-5,-5,-5,0.01), upper = c(5 setUp2 <- createBayesianSetup(likelihood2, lower = c(-5,-5,0.01), upper = c(5,5,30)) ``` -MCMC and marginal likelihood calculation +MCMC and marginal likelihood estimation ```{r, results = 'hide'} settings = list(iterations = 15000, message = FALSE) @@ -659,15 +705,15 @@ M2 ### Model comparison via Bayes factor -Bayes factor (need to reverse the log) +Bayes factor (need to invert the log) ```{r} exp(M1$ln.ML - M2$ln.ML) ``` -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. +BF \> 1 means that the evidence favors M1. See Kass, R. E. & Raftery, A. E. (1995) Bayes Factors. J. Am. Stat. Assoc. Amer Statist Assn, 90, 773-795. -Assuming equal prior weights on all models, we can calculate the posterior weight of M1 as +Assuming equal prior weights for all models, we can calculate the posterior weight of M1 as ```{r} exp(M1$ln.ML) / ( exp(M1$ln.ML) + exp(M2$ln.ML)) @@ -675,18 +721,16 @@ exp(M1$ln.ML) / ( exp(M1$ln.ML) + exp(M2$ln.ML)) If models have different model priors, multiply with the prior probabilities of each model. -### Model comparison via DIC - -The Deviance information criterion is a commonly applied method to summarize the fit of an MCMC chain. It can be obtained via +The Deviance Information Criterion is a commonly used method to summarize the fit of an MCMC chain. It can be calculated using ```{r} DIC(out1)$DIC DIC(out2)$DIC ``` -### Model comparison via WAIC +### Model Comparison via WAIC -The Watanabe--Akaike information criterion is another criterion for model comparison. To be able to calculate the WAIC, the model must implement a log-likelihood that density that allows to calculate the log-likelihood point-wise (the likelihood functions requires a "sum" argument that determines whether the summed log-likelihood should be returned). It can be obtained via +The Watanabe-Akaike Information Criterion is another criterion for model comparison. To compute the WAIC, the model must implement a log-likelihood density that allows the log-likelihood to be computed pointwise (the likelihood functions require a "sum" argument that determines whether the summed log-likelihood should be returned). It can be obtained via ```{r} # This will not work, since likelihood1 has no sum argument diff --git a/BayesianTools/vignettes/InterfacingAModel.Rmd b/BayesianTools/vignettes/InterfacingAModel.Rmd index 636490d..92f0fd4 100644 --- a/BayesianTools/vignettes/InterfacingAModel.Rmd +++ b/BayesianTools/vignettes/InterfacingAModel.Rmd @@ -7,7 +7,7 @@ vignette: > %\VignetteIndexEntry{Interfacing your model with R} \usepackage[utf8]{inputenc} %\VignetteEngine{knitr::rmarkdown} -abstract: "This tutorial discusses how to interface models written in other programming languages with R, so that they can be fit with BayesianTools" +abstract: "This tutorial discusses how to interface models written in other programming languages with R, so that they can be fit with BayesianTools" author: Florian Hartig editor_options: chunk_output_type: console @@ -25,10 +25,11 @@ set.seed(123) ## Step 1: Create a runModel(par) function -A basic requirement to allow calibration in BT is that we need to be able to execute the model with a given set of parameters. Strictly speaking, BT will not see the model as such, but requires a likelihood function with interface likelihood(par), where par is a vector, but in this function, you will then probably run the model with parameters par, where par could stay a vector, or be transformed to another format, e.g. data.frame, matrix or list. +To enable calibration in BT, it's essential to run the model with a specified set of parameters. In fact, BT does not see the model as such, but requires a likelihood function with the interface likelihood(par), where par is a vector, but in this function you will probably then run the model with the parameters par, where par could remain a vector or be transformed into some other format, e.g. data.frame, matrix or list. -What happens now depends on how your model is programmed - I have listed the steps in order of convenience / speed. If your model has never been interfaced with R you will likely have to move down to the last option. +What happens next depends on how your model is programmed. The following steps are arranged according to speed and convenience. If your model has not been interfaced with R before, you will most likely have to skip to the last option. +```{=tex} \begin{enumerate} \item Model in R, or R interface existing \item Model can be compiled and linked as a dll @@ -37,20 +38,20 @@ What happens now depends on how your model is programmed - I have listed the ste \item Model can be compiled as executable reads parameters via parameter file \item Model parameters are hard-coded in the executable \end{enumerate} - +``` ### Case 1 - model programmed in R -Usually, you will have to do nothing. Just make sure you can call your model a in +Typically, no action is required. Ensure that you are able to call your model. ```{r, eval = F} runMyModel(par) ``` -Typically this function will directly return the model outputs, so step 2 can be skipped. +Usually, this function returns model outputs directly, thus step 2 is unnecessary. ### Case 2 - compiled dll, parameters are set via dll interface -If you have your model prepared as a dll, or you can prepare it that way, you can use the \ref{https://stat.ethz.ch/R-manual/R-devel/library/base/html/dynload.html}{dyn.load()} function to link R to your model +If your model is already in the form of a DLL, or if you can prepare it as such, use the function \ref{https://stat.ethz.ch/R-manual/R-devel/library/base/html/dynload.html}{dyn.load()} to link R to your model. ```{r, eval = F} dyn.load(model) @@ -62,22 +63,21 @@ runMyModel(par){ } ``` -Again, if you implement this, you will also typically return the output directly via the dll and not write to file, which means that step 2 can be skipped. - -The tricky thing in this approach is that you have to code the interface to your dll, which technically means in most programming languages to set your variables as external or something similar, so that they can be accessed from the outside. How this works will depend on the programming language. +If you implement this process, you will also usually return the output directly via the `dll` and not write to a file, which means that step 2 can be skipped. +The problem with this approach is that you have to code the interface to your `dll`, which in most programming languages technically means setting your variables as external or something like that, so that they can be accessed from outside. The specifics of how this works will vary depending on the programming language being used. ### Case 3 - model programmed in C / C++, interfaced with RCPP -RCPP is a highly flexible environment to interface between R and C/C++. If your model is coded in C / C++, RCPP offers the most save and powerful way to connect with R (much more flexible than with command line or dll). +RCPP provides a highly flexible platform to interface between R and C/C++. RCPP provides a secure and powerful way to connect to R if your model is coded in C/C++ (much more flexible than using the command line or a `dll`). -However, doing the interface may need some adjustments to the code, and there can be technical problems that are difficult to solve for beginners. I do not recommend to attempt interfacing an existing C/C++ model unless you have RCPP experience, or at least very food C/C++ experience. +Nevertheless, code adjustments may be necessary for the interface, and beginners may find it challenging to resolve any technical issues. Attempting to interface an existing C/C++ model without prior experience using RCPP or at least substantial experience with C/C++ is not advisable. -Again, if you implement this, you will also typically return the output directly via the dll and not write to file, which means that step 2 can be skipped. +In addition, step 2 can be skipped if you are implementing this, as you will usually return the output directly via the `dll`, rather than writing to a file. ### Case 4 - compiled executable, parameters set via command line (std I/O) -If your model is written in a compiled or interpreted language, and accepts parameters via std I/O, wrapping is usually nothing more than writing the system call in an R function. An example would be +If your model is written in a compiled or interpreted language and accepts parameters via std I/O, wrapping is usually nothing more than writing the system call in an R function. For example ```{r, eval = F} runMyModel(par){ @@ -92,11 +92,11 @@ runMyModel(par){ } ``` -Note: If you have problems with the system command, try system2. If the model returns the output via std.out, you can catch this and convert it and skip step 2. If your model writes to file, go to step 2. +Note: If you encounter difficulties with the system command, try system2. If the model provides the output via std.out, you can catch and convert it and skip step 2. If your model writes the output to a file, proceed to step 2. ### Case 5 - compiled model, parameters set via parameter file or in any other method -Many models read parameters with a parameter file. In this case you want to do something like this +Many models use parameter files to read parameters. For this case, you want to do something like the following example ```{r, eval = F} runMyModel(par, returnData = NULL){ @@ -118,7 +118,7 @@ writeParameters(par){ } ``` -Depending on your problem, it can also make sense to define a setup function such as +For some problems, it may be useful to create a setup function, as in the example below. ```{r, eval = F} setUpModel <- function(parameterTemplate, site, localConditions){ @@ -130,31 +130,31 @@ setUpModel <- function(parameterTemplate, site, localConditions){ } ``` -How you do the write parameter function depends on the file format you use for the parameters. In general, you probably want to create a template parameter file that you use as a base and from which you change parameters +The way you write a parameter function depends on the file format you are using for the parameters. Usually, creating a template parameter file is recommended as a starting point, and then the parameters can be changed as required. -* If your parameter file is in an *.xml format*, check out the xml functions in R -* If your parameter file is in a *general text format*, the best option may be to create a template parameter file, place a unique string at the locations of the parameters that you want to replace, and then use string replace functions in R, e.g. [grep](https://stat.ethz.ch/R-manual/R-devel/library/base/html/grep.html) to replace this string. +- If your parameter file is in an *.xml format*, check out the xml functions in R. +- If your parameter file is in a *general text format*, the best option may be to create a template parameter file, place a unique string at the locations of the parameters that you want to replace, and then use string replace functions in R, e.g. [grep](https://stat.ethz.ch/R-manual/R-devel/library/base/html/grep.html) to replace this string. ### Case 6 - compiled model, parameters cannot be changed -You have to change your model code to achieve one of the former options. If the model is in C/C++, going directly to RCPP seems the best alternative. +To achieve one of the previous options, you must change your model code. If the model is in C/C++, the best alternative is to go directly to RCPP. -## Step 2: Reading back data +## Step 2: Read back data -If your model returns the output directly (which is highly preferable, ), you can skip this step. +If your model returns the output directly (which is highly preferable, ), you can skip this step. -For simple models, you might consider returning the model output directly with the runMyModel function. This is probably so for cases a) and b) above, i.e. model is already in R, or accepts parameters via command line. +You can skip this process if your model directly returns the output (which is highly recommended). -More complicated models, however, produce a large number of outputs and you typically don't need all of them. It is therefore more useful to make on or several separate readData or getDate function. The only two different cases I will consider here is +If you have a simple model, you might consider using the `runMyModel` function to return the model's output directly. This is suitable for cases in a) and b) namely, models that are already in R, or receive parameters via the command line. -* via dll / RCPP -* via file ouputs +In contrast, more complex models generate a large number of outputs. However, usually you do not need all of them. Therefore, it is more useful to create one or multiple separate `readData` or `getDate` functions. The only two cases I will consider here are -*Model is a dll* If the model is a dll file, the best thing would probably be to implement appropriate getData functions in the source code that can then be called from R. If your model is in C and in a dll, interfacing this via RCPP would probably be easier, because you can directly return R dataframes and other data structures. +- via dll / RCPP +- via file ouputs +*Model is a dll* If the model is a `dll` file, it would probably be best to implement appropriate `getData` functions in the source code that can then be called from R. If your model is in C and in a `dll`, interfacing it via RCPP would probably be easier because you can return R dataframes and other data structures directly. -*Model writes file output* If the model writes file output, write a getData function that reads in the model outputs and returns the data in the desired format, typically the same that you would use to represent your field data. - +*Model writes file output* If the model generates file output, create a `getData` function that reads the model outputs and returns the data in the desired format, typically the same format you would use to represent your field data. ```{r, eval = F} getData(type = X){ @@ -167,12 +167,9 @@ getData(type = X){ } ``` +\subsection{Testing the approach} -\subsection{Testing the approach} - - -From R, you should now be able to do something like that - +You should now see an example from R that looks like this ```{r, eval = F} par = c(1,2,3,4 ..) @@ -184,55 +181,51 @@ output <- getData(type = DesiredType) plot(output) ``` - - - ## Step 3 (optional) - creating an R package from your code -The last step is optional, but we recommend that you take it from the start, because there is really no downside to it. You work with R code in several files that you run by hand, or diretly put all code into an R package. Creating and managing R packages is very easy, and it's easier to pass on your code because everything, including help, is in on package. To create an R package, follow the tutorial \href{http://biometry.github.io/APES/R/R70-PackageDevelopment.html}{here}. Remember to create a good documentation using Roxygen. - +Although the final step is optional, we recommend doing it from the beginning because there is really no downside. You can work with R code in multiple files that can be run manually or incorporated into an R package directly. Creating and managing R packages is a straightforward process and makes it simpler to share your code, as everything, including help guides, is in one package. To create an R package, please refer to the tutorial \href{http://biometry.github.io/APES/R/R70-PackageDevelopment.html}{here}. Please remember to write good documentation using Roxygen. # Speed optimization and parallelization -For running sensitivity analyses or calibrations, runtime is often an issue. Before you parallelize, make sure your model is as fast as possible. +Runtime is often a concern when performing sensitivity analyses or calibrations. Ensure that your model has been optimized for maximum speed before parallelization. ## Easy things -* Are you compiling with maximum optimization (e.g. -o3 in cpp) -* If you have a spin-up phase, could you increase the time-step during this phase? -* Could you increase the time step generally -* Do you write unnecessary outputs that you could turn off (harddisk I/O is often slow)? +- Are you compiling with maximum optimization (e.g. -o3 in cpp) +- If there is a spin-up phase, consider increasing the time step during that phase. +- Consider increasing the time step in general. +- Are you producing unnecessary outputs that could be turned off to reduce hard disk I/O, which is often slow? ## Difficult things -* Make the model directly callable (RCPPor dll) to avoid harddisk I/O -* Is it possible to reduce initialization time (not only spin-up, but also for reading in the forcings / drivers) by avoid ending the model executable after each run, but rather keep it "waiting" for a new run. -* Code optimization: did you use a profiler? Read up on code optimization -* Check for unnecessary calculations in your code / introduce compiler flags where appropriate +- Make the model directly callable (RCPP or dll) to avoid harddisk I/O +- Is it possible to reduce the initialization time (not only for spin-up, but also for reading the forcings / drivers) by not terminating the model executable after each run, but keeping it "waiting" for a new run? +- Code optimization: did you use a profiler? Read up on code optimization +- Check for unnecessary calculations in your code / introduce compiler flags where appropriate ## Parallelization -A possibility to speed up the run time of your model is to run it on multiple cores (CPU's). To do so, you have two choices: +One way to speed up the runtime of your model is to run it on multiple cores (CPUs). To do so, you have two choices: -1. Parallelize the model itself -2. Parallelize the model call, so that BT can do several model evaluations in parallel +1. Parallelize the model itself +2. Parallelize the model call so that BT can perform multiple model evaluations in parallel. -Which of the two makes more sense depends a lot on your problem. To parallelize the model itself will be interesting in particular for very large models, which could otherwise not be calibrated with MCMCs. However, this approach will typically require to write parallel C/C++ code and require advanced programming skills, which is the reason why we will not further discuss it here. +Which of the two makes more sense depends a lot on your problem. Parallelizing the model itself will be especially interesting for very large models that could not otherwise be calibrated with MCMCs. However, this approach typically requires writing parallel C/C++ code and advanced programming skills, which is why we will not discuss it further here. -The usual advice in parallel computing is anyway to parallelize the outer loops first, to minimize communication overhead, which would suggest to start with parallelizing model evaluations. This is also much easier to program. Even within this, there are two levels of parallelization possible: +The usual advice in parallel computing anyway is to parallelize the outer loops first to minimize the communication overhead, which would suggest starting with parallelizing the model evaluations. This is also much easier to program. Even within that, there are two levels of parallelization possible: -1. Parallelize the call of several MCMC / SMC samplers -2. Parallelize within the MCMC / SMC samplers +1. Parallelize the call of multiple MCMC / SMC samplers. +2. Parallelize within the MCMC / SMC samplers -Currently, BT only supports parallelization within MCMCs / SMCs, but it easy to also implement between sampler parallelization by hand. Both approaches are describe below. +Currently, BT only supports parallelization within MCMCs / SMCs, but it easy to also implement between sampler parallelization by hand. Both approaches are describe below. -### Within sampler parallelization +### Within sampler parallelization -Within-sampler parallelization is particular useful for algorithms that can use a large number of cores in parallel, e.g. sensitivity analyses or SMC sampling. For the MCMCs, it depends on the settings and the algorithm how much parallelization they can make use of. In general, MCMCs are Markovian, as the name says, i.e. the set up a chain of sequential model evaluations, and those calls can therefore not be fully parallelized. However, a number of MCMCs in the BT package uses MCMC algorithms that can be partly parallelized, in particular the population MCMC algorithms DE/DEzs/DREAM/DREAMzs. For all these cases, BT will automatically use parallelization of the BT setup indicates that this is implemented. +Within-sampler parallelization is particularly useful for algorithms that can use a large number of cores in parallel, such as sensitivity analysis or SMC sampling. For MCMCs, the amount of parallelization that can be used depends on the settings and the algorithm. In general, MCMCs are, as the name implies, Markovian, i.e., they set up a chain of sequential model evaluations, and these calls cannot be fully parallelized. However, a number of MCMCs in the BT package uses MCMC algorithms that can be partly parallelized, in particular the population MCMC algorithms DE/DEzs/DREAM/DREAMzs. In all these cases, BT will automatically use parallelization of the BT setup to indicate that it is implemented. -How to do this? A first requirement to do so is to to have your model wrapped into an R function (see PREVIOUS SECTION). If that is the case, R offers a number of options to run functions in parallel. The easiest is to use the parallel package that comes with the R core. For other packages, see the internet and the CRAN task view on [High Performance Computing](https://CRAN.R-project.org/view=HighPerformanceComputing) +How to do this? A first requirement is that your model is wrapped in an R function (see PREVIOUS SECTION). Once that is the case, R provides a number of ways to run functions in parallel. The easiest is to use the parallel package that comes with the R core. Other packages can be found on the Internet and in the CRAN task view at [High Performance Computing](https://CRAN.R-project.org/view=HighPerformanceComputing) -As an example, assume we have the following, very simple model: +As an example, suppose we have the following very simple model: ```{r} mymodel<-function(x){ @@ -241,7 +234,7 @@ mymodel<-function(x){ } ``` -To start a parallel computation, we will first need to create a cluster object. Here we will initiate a cluster with 2 CPU's. +To start a parallel computation, we must first create a cluster object. Here we start a cluster with 2 CPUs. ```{r, eval = F} @@ -255,18 +248,18 @@ runParallel<- function(parList){ runParallel(c(1,2)) ``` -You could use this principle to build your own parallelized likelihood. However, something very similar to the previous loop is automatized in BayesianTools. You can directly create a parallel model evaluation function with the function generateParallelExecuter, or alternatively directly in the createBayesianSetup +You could use this principle to build your own parallelized likelihood. However, something very similar to the previous loop is automated in `BayesianTools`. You can create a parallel model evaluation function directly with the `generateParallelExecuter` function, or alternatively directly in the `createBayesianSetup` function. ```{r, eval = F} library(BayesianTools) parModel <- generateParallelExecuter(mymodel) ``` -If your model is tread-safe, you should be able to run this out of the box. I therefore recommend using the hand-coded paraellelization only of the model is not thread-safe. +If your model is thread-safe, you should be able to run it out of the box. Therefore, I recommend using the hand-coded parallelization only if your model is not thread-safe. ### Running several MCMCs in parallel -Additionally to the within-chain parallelization, you can also run several MCMCs in parallel, and combine them later to a single McmcSamplerList +In addition to within-chain parallelization, you can also run multiple MCMCs in parallel and later combine them into a single `McmcSamplerList`. ```{r} library(BayesianTools) @@ -284,17 +277,8 @@ res <- createMcmcSamplerList(list(out1, out2)) plot(res) ``` - - ### Thread safety -Thread safety quite generally means that you can execute multiple instances of your code on your hardware. There are various things that can limit Thread safety, for example - -* writing outputs to file (several threads might write to the same file at the same time) - - - - - - +Thread safety generally means that you can run multiple instances of your code on your hardware. There are several things that can limit thread safety, such as +- Writing output to a file (multiple threads can write to the same file at the same time)