MCMCEnsemble: MCMC ensemble sampler

In SandaD/MCMCEnsembleSampler: MCMC sampler

Description

Ensemble Markov Chain Monte Carlo sampler with different strategies to generate proposals. Either the stretch move as proposed by Goodman and Weare (2010), or a differential evolution jump move similar to Braak and Vrugt (2008).

Usage

MCMCEnsemble(f, lower.inits, upper.inits, max.iter, n.walkers = 10 *
  length(lower.inits), method = c("stretch", "differential.evolution"),
  coda = FALSE, ...)

Arguments

f	function that returns a value proportional to the log probability density to sample from.
lower.inits	vector specifying for each parameters the lower value the initial distibution
upper.inits	vector specifying for each parameters the upper value the initial distibution
max.iter	maximum number of function evaluations
n.walkers	number of walkers (ensemble size)
method	method for proposal generation, either "stretch", or "differential.evolution". The first letter is sufficient.
coda	logical. Should the samples be returned as link[coda]{mcmc.list} object?
...	further arguments passed to f

Value

if coda==FALSE a list with

• samples A three dimensional array of samples with dimensions walker x generation x parameter

• log.p A matrix with the log density evaluate for aeach walker at each generation.

if coda==TRUE a list with

• samples A object of class link[coda]{mcmc.list} containing all samples.

• log.p A matrix with the log density evaluate for aeach walker at each generation.

References

Goodman, J. and Weare, J. (2010) Ensemble samplers with affine invariance. Communications in Applied Mathematics and Computational Science, 5(1), 65-80.

Braak, C. J. F. ter and Vrugt, J. A. (2008) Differential Evolution Markov Chain with snooker updater and fewer chains. Statistics and Computing, 18(4), 435-446.

Examples

## a log-pdf to sample from
p.log <- function(x) {
    B <- 0.03                              # controls 'bananacity'
    -x[1]^2/200 - 1/2*(x[2]+B*x[1]^2-100*B)^2
}

## use stretch move
res1 <- MCMCEnsemble(p.log, lower.inits=c(a=0, b=0), upper.inits=c(a=1, b=1),
                     max.iter=3000, n.walkers=10, method="s")
str(res1)


## use stretch move, return samples as 'coda' object
res2 <- MCMCEnsemble(p.log, lower.inits=c(a=0, b=0), upper.inits=c(a=1, b=1),
                     max.iter=3000, n.walkers=10, method="s", coda=TRUE)

summary(res2$samples)
plot(res2$samples)


## use different evolution move, return samples as 'coda' object
res3 <- MCMCEnsemble(p.log, lower.inits=c(a=0, b=0), upper.inits=c(a=1, b=1),
                     max.iter=3000, n.walkers=10, method="d", coda=TRUE)

summary(res3$samples)
plot(res3$samples)

SandaD/MCMCEnsembleSampler documentation built on May 9, 2019, 12:25 p.m.