MCMC {adaptMCMC} | R Documentation |
(Adaptive) Metropolis Sampler
Description
Implementation of the robust adaptive Metropolis sampler of Vihola (2012).
Usage
MCMC(p, n, init, scale = rep(1, length(init)),
adapt = !is.null(acc.rate), acc.rate = NULL, gamma = 2/3,
list = TRUE, showProgressBar=interactive(), n.start = 0, ...)
Arguments
p |
function that returns a value proportional to the log probability
density to sample from. Alternatively it can be a function that returns a list
with at least one element named |
n |
number of samples. |
init |
vector with initial values. |
scale |
vector with the variances or covariance matrix of the jump distribution. |
adapt |
if |
acc.rate |
desired acceptance rate (ignored if |
gamma |
controls the speed of adaption. Should be between 0.5 and 1. A lower gamma leads to faster adaption. |
list |
logical. If |
showProgressBar |
logical. If |
n.start |
iteration where the adaption starts. Only internally used. |
... |
further arguments passed to |
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.
In some cases the function p
may not only calculate the log
density but return a list containing also other values. For example
if p
is a log posterior one may be also interested to store
the corresponding prior and likelihood values. The function must
either return always a scalar or always a list, however, the length of
the list may vary.
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 |
log.p |
vector with the (unnormalized) log density for each sample |
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 |
.
extra.values |
A list containing additional return values provided by
|
Note
Due to numerical errors it may happen that the computed
covariance matrix is not positive definite. In such a case the nearest positive
definite matrix is calculated with nearPD()
from the package Matrix.
Author(s)
Andreas Scheidegger, andreas.scheidegger@eawag.ch or scheidegger.a@gmail.com.
Thanks to David Pleydell, Venelin, and Umberto Picchini for spotting
errors and providing improvements. Ian Taylor implemented the usage of
adapt_S
which lead to a nice speedup.
References
Vihola, M. (2012) Robust adaptive Metropolis algorithm with coerced acceptance rate. Statistics and Computing, 22(5), 997-1008. doi:10.1007/s11222-011-9269-5.
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 last 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)
## ----------------------
## function returning extra information in a list
p.log.list <- function(x) {
B <- 0.03 # controls 'bananacity'
log.density <- -x[1]^2/200 - 1/2*(x[2]+B*x[1]^2-100*B)^2
result <- list(log.density=log.density)
## under some conditions one may want to return other infos
if(x[1]<0) {
result$message <- "Attention x[1] is negative!"
result$x <- x[1]
}
result
}
samp.list <- MCMC(p.log.list, n=200, init=c(0, 1), scale=c(1, 0.1),
adapt=TRUE, acc.rate=0.234)
## the additional values are stored under `extras.values`
head(samp.list$extras.values)