adaptMCMC (version 1.0.2)

MCMC: (Adaptive) Metropolis Sampler

Description

Implementation of the robust adaptive Metropolis sampler of Vihola (2011).

Usage

MCMC(p, n, init, scale = rep(1, length(init)), log = TRUE,
    adapt = !is.null(acc.rate), acc.rate = NULL, gamma = 0.55,
    list = TRUE, n.start = 0, ...)

Arguments

p

function that returns the (log) probability density to sample from. Must have two or more dimensions.

n

number of samples.

init

vector with initial values.

scale

vector with the variances or covariance matrix of the jump distribution.

log

logical. If TRUE, a log density is expected from p (strongly recommended).

adapt

if TRUE, adaptive sampling is used, if FALSE classic metropolis sampling, if a positive integer the adaption stops after adapt iterations.

acc.rate

desired acceptance rate (ignored if adapt=FALSE)

gamma

controls the speed of adaption. Should be between 0.5 and 1. A lower gamma leads to faster adaption.

list

logical. If TRUE a list is returned otherwise only a matrix with the samples.

n.start

iteration where the adaption starts. Only internally used.

...

further arguments passed to p.

Value

• If list=FALSE a matrix is with the samples. If list=TRUE a list is returned with the following components:
• samples
  matrix with samples
• n.sample
  number of generated samples
• acceptance.rate
  acceptance rate
• adaption
  Either logical if adaption was used or not, or the number of adaption steps.
• sampling.parameters
  a list with further sampling parameters. Mainly used by MCMC.add.samples().

Details

The algorithm tunes the covariance matrix of the (normal) jump distribution to achieve the desired acceptance rate. Classic (non-adaptive) Metropolis sampling can be obtained by setting adapt=FALSE. Note, due to the calculation for the adaption steps the sampler is rather slow. However, with a suitable jump distribution good mixing can be observed with less samples. This is crucial if the computation of p is slow. The current implementation does not work for one-dimensional distributions.

References

Vihola, M. (2011) Robust adaptive Metropolis algorithm with coerced acceptance rate. Statistics and Computing. [online] http://www.springerlink.com/content/672270222w79h431/, (Accessed December 8, 2011)

See Also

MCMC.parallel, MCMC.add.samples

Examples

## ----------------------
## Banana shaped distribution

## 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
}


## ----------------------
## generate samples

## 1) non-adaptive sampling
samp.1 <- MCMC(p.log, n=200, init=c(0, 1), scale=c(1, 0.1),
               adapt=FALSE)

## 2) adaptive sampling
samp.2 <- MCMC(p.log, n=200, init=c(0, 1), scale=c(1, 0.1),
               adapt=TRUE, acc.rate=0.234)


## ----------------------
## summarize results

str(samp.2)
summary(samp.2$samples)

## covariance of jump distribution
samp.2$cov.jump


## ----------------------
## plot density and samples

x1 <- seq(-15, 15, length=80)
x2 <- seq(-15, 15, length=80)
d.banana <- matrix(apply(expand.grid(x1, x2), 1,  p.log), nrow=80)

par(mfrow=c(1,2))
image(x1, x2, exp(d.banana), col=cm.colors(60), asp=1, main="no adaption")
contour(x1, x2, exp(d.banana), add=TRUE, col=gray(0.6))
lines(samp.1$samples, type='b', pch=3)

image(x1, x2, exp(d.banana), col=cm.colors(60), asp=1, main="with adaption")
contour(x1, x2, exp(d.banana), add=TRUE, col=gray(0.6))
lines(samp.2$samples, type='b', pch=3)