metrop package:mcmc R Documentation
_M_e_t_r_o_p_o_l_i_s _A_l_g_o_r_i_t_h_m
_D_e_s_c_r_i_p_t_i_o_n:
Markov chain Monte Carlo for continuous random vector using a
Metropolis algorithm.
_U_s_a_g_e:
metrop(obj, initial, nbatch, blen = 1, nspac = 1, scale = 1, outfun,
debug = FALSE, ...)
_A_r_g_u_m_e_n_t_s:
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’.
_D_e_t_a_i_l_s:
Runs a “random-walk” Metropolis algorithm 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, ...) > 0’. Description of a complete MCMC
analysis (Bayesian logistic regression) using this function can be
found in the vignette ‘demo’ ().
_V_a_l_u_e:
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’ ().
_W_a_r_n_i_n_g:
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.
_E_x_a_m_p_l_e_s:
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])