Parallelization of the likelihood evaluations

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.

Algorithms in the BayesianTools package can make use of parallel computing if this option is specified in the BayesianSetup. Note that currently, parallelization is used by the following algorithms: SMC, DEzs and DREAMzs sampler. It can also be used through the BayesianSetup with the functions of the sensitivity package.

Here are some more details about the parallelization.

## Definition of likelihood function
likelihood <- function(matrix){
    # Calculate likelihood in parallel
    # Return vector of likelihood values
}

## Create Bayesian Setup
BS <- createBayesianSetup(likelihood, parallel = "external", ...)

## Run MCMC
runMCMC(BS, sampler = "SMC", ...)

## n = Number of cores
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)

## give runMCMC a matrix with n rows of proposals as startValues or sample n times from the previous created sampler
out <- runMCMC(bayesianSetup, settings = list(iterations = 1000))

### Create cluster with n cores
cl <- parallel::makeCluster(n)

## Definition of the likelihood
likelihood  <- function(X) sum(dnorm(c(1:10), mean = X, log = T))

## Definition of the likelihood which will be calculated in parallel. Instead of the parApply function, we could also define a costly parallelized likelihood
pLikelihood <- function(param) parallel::parApply(cl = cl, X = param, MARGIN = 1, FUN = likelihood)

## export functions, dlls, libraries
# parallel::clusterEvalQ(cl, library(BayesianTools))
parallel::clusterExport(cl, varlist = list(likelihood))

## create BayesianSetup
bayesianSetup <- createBayesianSetup(pLikelihood, lower = -10, upper = 10, parallel = 'external')

## For this case we want to parallelize the internal chains, therefore we create a n row matrix with startValues, if you parallelize a model in the likelihood, do not set a n*row Matrix for startValue
settings = list(iterations = 100, nrChains = 1, startValue = bayesianSetup$prior$sampler(n))

## runMCMC
out <- runMCMC(bayesianSetup, settings, sampler = "DEzs")

### Create cluster with n cores
cl <- parallel::makeCluster(n)

## export your model
# parallel::clusterExport(cl, varlist = list(complexModel))

## Definition of the likelihood
likelihood  <- function(param) {
  # ll <- complexModel(param)
  # return(ll)
} 

## create BayesianSetup and settings
bayesianSetup <- createBayesianSetup(likelihood, lower = -10, upper = 10, parallel = 'external')
settings = list(iterations = 100, nrChains = 1)

## runMCMC
out <- runMCMC(bayesianSetup, settings)

### Definition of likelihood function
x <- c(1:10)
likelihood <- function(param) return(sum(dnorm(x, mean = param, log = T)))

## Create BayesianSetup and settings
bayesianSetup <- createBayesianSetup(likelihood, lower = -10, upper = 10, parallel = F)
settings = list(iterations = 100000)

## Start cluster with n cores for n chains and export BayesianTools library
cl <- parallel::makeCluster(n)
parallel::clusterEvalQ(cl, library(BayesianTools))

## calculate parallel n chains, for each chain the likelihood will be calculated on one core
out <- parallel::parLapply(cl, 1:n, fun = function(X, bayesianSetup, settings) runMCMC(bayesianSetup, settings, sampler = "DEzs"), bayesianSetup, settings)

## Combine the chains
out <- createMcmcSamplerList(out)

Creating priors

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-define prior, see ‘?createPrior’
• Create a prior from a previous MCMC sample

The prior we choose depends on the prior information we have. For example, if we have no prior information, we can choose a uniform prior. The normal distribution is often used to model a wide range of phenomena in statistics, such as the distribution of heights or weights in a population. Beta distribution, on the other hand, is defined on the interval 0, 1. It is often used to model random variables that represent proportions, probabilities or other values that are constrained to lie within this interval.

$\color{darkgreen}{\text{createPrior}}$	$\color{darkgreen}{\text{createBetaPrior}}$	$\color{darkgreen}{\text{createUniformPrior}}$	$\color{darkgreen}{\text{createTruncatedNormalPrior}}$
Any density provided by the user	Beta density	Uniform density	Normal density
createPrior(density, sampler, lower, upper, best)	createBetaPrior(a, b, lower, upper)	createUniformPrior(lower, upper, best)	createTruncatedNormalPrior(mean, sd, lower, upper).

# Create a BayesianSetup
ll <- generateTestDensityMultiNormal(sigma = "no correlation")
bayesianSetup = createBayesianSetup(likelihood = ll, lower = rep(-10, 3), upper = rep(10, 3))

settings = list(iterations = 2500,  message = FALSE)
out <- runMCMC(bayesianSetup = bayesianSetup, settings = settings)


newPrior = createPriorDensity(out, method = "multivariate", eps = 1e-10, lower = rep(-10, 3), upper =  rep(10, 3), best = NULL)
bayesianSetup <- createBayesianSetup(likelihood = ll, prior = newPrior)

settings = list(iterations = 1000,  message = FALSE)
out <- runMCMC(bayesianSetup = bayesianSetup, settings = settings)

The 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 specifying sampler = "Metropolis" in the runMCMC().

The subsequent code provides an overview of the default settings of the MH sampler.

The following code gives an overview about the default settings of the MH sampler.

applySettingsDefault(sampler = "Metropolis")

## $sampler
## [1] "Metropolis"
## 
## $startValue
## NULL
## 
## $iterations
## [1] 10000
## 
## $burnin
## [1] 0
## 
## $thin
## [1] 1
## 
## $consoleUpdates
## [1] 100
## 
## $parallel
## NULL
## 
## $message
## [1] TRUE
## 
## $optimize
## [1] TRUE
## 
## $proposalGenerator
## NULL
## 
## $adapt
## [1] FALSE
## 
## $adaptationInterval
## [1] 500
## 
## $adaptationNotBefore
## [1] 3000
## 
## $DRlevels
## [1] 1
## 
## $proposalScaling
## NULL
## 
## $adaptationDepth
## NULL
## 
## $temperingFunction
## NULL
## 
## $gibbsProbabilities
## NULL
## 
## $nrChains
## [1] 1
## 
## $runtime
## [1] 0
## 
## $sessionInfo
## R version 4.4.2 (2024-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.1 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] BayesianTools_0.1.8 rmarkdown_2.29     
## 
## loaded via a namespace (and not attached):
##  [1] Matrix_1.7-1         jsonlite_1.8.9       compiler_4.4.2      
##  [4] Rcpp_1.0.13-1        stringr_1.5.1        DHARMa_0.4.7        
##  [7] jquerylib_0.1.4      splines_4.4.2        boot_1.3-31         
## [10] yaml_2.3.10          fastmap_1.2.0        lattice_0.22-6      
## [13] coda_0.19-4.1        R6_2.5.1             Brobdingnag_1.2-9   
## [16] knitr_1.49           MASS_7.3-61          tmvtnorm_1.6        
## [19] nloptr_2.1.1         maketools_1.3.1      minqa_1.2.8         
## [22] bslib_0.8.0          rlang_1.1.4          cachem_1.1.0        
## [25] stringi_1.8.4        xfun_0.49            IDPmisc_1.1.21      
## [28] sass_0.4.9           sys_3.4.3            cli_3.6.3           
## [31] magrittr_2.0.3       digest_0.6.37        grid_4.4.2          
## [34] gmm_1.8              mvtnorm_1.3-2        sandwich_3.1-1      
## [37] lme4_1.1-35.5        lifecycle_1.0.4      nlme_3.1-166        
## [40] evaluate_1.0.1       glue_1.8.0           numDeriv_2016.8-1.1 
## [43] zoo_1.8-12           buildtools_1.0.0     bridgesampling_1.1-2
## [46] stats4_4.4.2         tools_4.4.2          htmltools_0.5.8.1

Here are some examples of how to apply different settings. Activate individual options or combinations as demonstrated.

Standard MH 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.

Metropolis, N., A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller (1953). Equation of state calculations by fast computing machines. The journal of chemical physics 21 (6), 1087 - 1092.

settings <- list(iterations = iter, adapt = F, DRlevels = 1, gibbsProbabilities = NULL, temperingFunction = NULL, optimize = F, message = FALSE)
out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings)
plot(out) 

settings <- list(iterations = iter, adapt = F, DRlevels = 1, gibbsProbabilities = NULL, temperingFunction = NULL, optimize = T, message = FALSE)
out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings)
plot(out) 

settings <- list(iterations = iter, adapt = T, DRlevels = 1, gibbsProbabilities = NULL, temperingFunction = NULL, optimize = T,  message = FALSE)
out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings)
plot(out) 

settings <- list(iterations = iter, adapt = F, DRlevels = 2, gibbsProbabilities = NULL, temperingFunction = NULL, optimize = T,  message = FALSE)
out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings)
plot(out) 

settings <- list(iterations = iter, adapt = T, DRlevels = 2, gibbsProbabilities = NULL, temperingFunction = NULL, optimize = T,  message = FALSE)
out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings)
plot(out) 

settings <- list(iterations = iter, adapt = T, DRlevels = 1, gibbsProbabilities = c(1,0.5,0), temperingFunction = NULL, optimize = T,  message = FALSE)
out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings)
plot(out) 

temperingFunction <- function(x) 5 * exp(-0.01*x) + 1
settings <- list(iterations = iter, adapt = F, DRlevels = 1, gibbsProbabilities = c(1,1,0), temperingFunction = temperingFunction, optimize = T,  message = FALSE)
out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings)
plot(out) 

Differential Evolution MCMC

The BT package implements two differential evolution MCMCs. When in doubt, use the DEzs option.

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 samplers 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.

settings <- list(iterations = iter,  message = FALSE)
out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "DE", settings = settings)
plot(out) 

The second is the Differential Evolution MCMC with snooker update and sampling from past states, 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.

settings <- list(iterations = iter,  message = FALSE)
out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "DEzs", settings = settings)
plot(out) 

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.

This sampler is largely based on the DE sampler with some significant changes:

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.

settings <- list(iterations = iter,  message = FALSE)
out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "DREAM", settings = settings)
plot(out) 

The second implementation uses the same extension as the DEzs sampler. Namely sampling from past states and a snooker update. Again, this extension allows the use of fewer chains and parallel computing.

Again, if in doubt, you should prefer “DREAMzs”.

settings <- list(iterations = iter,  message = FALSE)
out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "DREAMzs", settings = settings)
plot(out) 

T-walk

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.

settings = list(iterations = iter,  message = FALSE) 

out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Twalk", settings = settings)

Rejection samling

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.

settings <- list(initialParticles = iter, iterations= 1)
out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "SMC", settings = settings)
plot(out) 

Sequential Monte Carlo (SMC)

The more sophisticated option is to use the implemented SMC, which is basically a particle filter that applies multiple filter steps.

settings <- list(initialParticles = iter, iterations= 10)
out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "SMC", settings = settings)
plot(out) 

Note that using a number for initialParticles requires that the bayesianSetup includes the option to sample from the prior.

Model comparison via Bayes factor

Bayes factor (need to invert the log)

exp(M1$ln.ML - M2$ln.ML)

## [1] 2.613552e+21

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 for all models, we can calculate the posterior weight of M1 as

exp(M1$ln.ML) / ( exp(M1$ln.ML) + exp(M2$ln.ML))

## [1] 1

If models have different model priors, multiply with the prior probabilities of each model.

The Deviance Information Criterion is a commonly used method to summarize the fit of an MCMC chain. It can be calculated using

DIC(out1)$DIC

## [1] 269.7765

DIC(out2)$DIC

## [1] 364.9808

Model Comparison via WAIC

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

# This will not work, since likelihood1 has no sum argument
# WAIC(out1)

# likelihood with sum argument
likelihood3 <- function(param, sum = TRUE){
    pred <- param[1] + param[2]*x + param[3] * x^2
    singlelikelihoods <- dnorm(y, mean = pred, sd = 1/(param[4]^2), log = T)
    return(if (sum == TRUE) sum(singlelikelihoods) else singlelikelihoods)  
}
setUp3 <- createBayesianSetup(likelihood3, lower = c(-5,-5,-5,0.01), upper = c(5,5,5,30))
out3 <- runMCMC(bayesianSetup = setUp3, sampler = "Metropolis", settings = settings)

WAIC(out3)

## $WAIC1
## [1] 1104035854
## 
## $WAIC2
## [1] 1104035854
## 
## $lppd
## [1] -552017927
## 
## $pWAIC1
## [1] 1.848716e-05
## 
## $pWAIC2
## [1] 1.851485e-05