R: Metropolis Algorithm

metrop {mcmc}

R Documentation

Metropolis Algorithm

Description

Markov chain Monte Carlo for continuous random vector using a Metropolis algorithm.

Usage

metrop(obj, initial, nbatch, blen = 1, nspac = 1, scale = 1, outfun,
    debug = FALSE, ...)

Arguments

`obj`	an R function that evaluates the log unnormalized probability density of the desired equilibrium distribution of the Markov chain. First argument is the state vector of the Markov chain. Other arguments arbitrary and taken from the `...` arguments of this function. Should return `- Inf` for points of the state space having probability zero under the desired equilibrium distribution. Alternatively, an object of class `"metropolis"` from a previous run can be supplied, in which case any missing arguments (including the log unnormalized density function) are taken from this object (up until version 0.7-2 this was incorrect with respect to the `debug` argument, now it applies to it too).
`initial`	a real vector, the initial state of the Markov chain.
`nbatch`	the number of batches.
`blen`	the length of batches.
`nspac`	the spacing of iterations that contribute to batches.
`scale`	controls the proposal step size. If scalar or vector, the proposal is `x + scale * z` where `x` is the current state and `z` is a standard normal random vector. If matrix, the proposal is `x + scale %*% z`.
`outfun`	controls the output. If a function, then the batch means of `outfun(state, ...)` are returned. If a numeric or logical vector, then the batch means of `state[outfun]` (if this makes sense). If missing, the the batch means of `state` are returned.
`debug`	if `TRUE` extra output useful for testing.
`...`	additional arguments for `obj` or `outfun`.

Details

Runs a “random-walk” Metropolis algorithm, terminology introduced by Tierney (1994), with multivariate normal proposal producing a Markov chain with equilibrium distribution having a specified unnormalized density. Distribution must be continuous. Support of the distribution is the support of the density specified by argument obj. The initial state must satisfy obj(state, ...) > - Inf. Description of a complete MCMC analysis (Bayesian logistic regression) using this function can be found in the vignette demo (../doc/demo.pdf).

Suppose the function coded by the log unnormalized function (either obj or obj$lud) is actually a log unnormalized density, that is, if w denotes that function, then exp(w) integrates to some value strictly between zero and infinity. Then the metrop function always simulates a reversible, Harris ergodic Markov chain having the equilibrium distribution with this log unnormalized density. The chain is not guaranteed to be geometrically ergodic. In fact it cannot be geometrically ergodic if the tails of the log unnormalized density are sufficiently heavy. The morph.metrop function deals with this situation.

Value

an object of class "mcmc", subclass "metropolis", which is a list containing at least the following components:

`accept`	fraction of Metropolis proposals accepted.
`batch`	`nbatch` by `p` matrix, the batch means, where `p` is the dimension of the result of `outfun` if `outfun` is a function, otherwise the dimension of `state[outfun]` if that makes sense, and the dimension of `state` when `outfun` is missing.
`initial`	value of argument `initial`.
`final`	final state of Markov chain.
`initial.seed`	value of `.Random.seed` before the run.
`final.seed`	value of `.Random.seed` after the run.
`time`	running time of Markov chain from `system.time()`.
`lud`	the function used to calculate log unnormalized density, either `obj` or `obj$lud` from a previous run.
`nbatch`	the argument `nbatch` or `obj$nbatch`.
`blen`	the argument `blen` or `obj$blen`.
`nspac`	the argument `nspac` or `obj$nspac`.
`outfun`	the argument `outfun` or `obj$outfun`.

Description of additional output when debug = TRUE can be found in the vignette debug (../doc/debug.pdf).

Warning

If outfun is missing or not a function, then the log unnormalized density can be defined without a ... argument and that works fine. One can define it starting ludfun <- function(state) and that works or ludfun <- function(state, foo, bar), where foo and bar are supplied as additional arguments to metrop.

If outfun is a function, then both it and the log unnormalized density function can be defined without ... arguments if they have exactly the same arguments list and that works fine. Otherwise it doesn't work. Start the definitions ludfun <- function(state, foo) and outfun <- function(state, bar) and you get an error about unused arguments. Instead start the definitions ludfun <- function(state, foo, ...) and outfun <- function(state, bar, ...), supply foo and bar as additional arguments to metrop, and that works fine.

In short, the log unnormalized density function and outfun need to have ... in their arguments list to be safe. Sometimes it works when ... is left out and sometimes it doesn't.

Of course, one can avoid this whole issue by always defining the log unnormalized density function and outfun to have only one argument state and use global variables (objects in the R global environment) to specify any other information these functions need to use. That too follows the R way. But some people consider that bad programming practice.

References

Tierney, L. (1994) Markov chains for exploring posterior distributions (with discussion). Annals of Statistics 22 1701–1762.

Examples

h <- function(x) if (all(x >= 0) && sum(x) <= 1) return(1) else return(-Inf)
out <- metrop(h, rep(0, 5), 1000)
out$accept
# acceptance rate too low
out <- metrop(out, scale = 0.1)
out$accept
# acceptance rate o. k. (about 25 percent)
plot(out$batch[ , 1])
# but run length too short (few excursions from end to end of range)
out <- metrop(out, nbatch = 1e4)
out$accept
plot(out$batch[ , 1])
hist(out$batch[ , 1])

[Package mcmc version 0.9-4 Index]