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 Statistics22 1701–1762.
See Also
morph.metrop
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])
providing 
.