metrop(obj, initial, nbatch, blen = 1, nspac = 1, scale = 1, outfun,
debug = FALSE, ...)x + scale * z where x is
the current state and z is a standard normal random vector.
If matrix, the proposal is outfun(state, ...) are returned. If a numeric
or logical vector, then the batch means of state[outfun]
(if this makes sense). If missingTRUE extra output useful for testing.obj or outfun."mcmc", subclass "metropolis",
which is a list containing at least the following components: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..Random.seed before the run..Random.seed after the run.system.time().obj or obj$lud from a previous run.nbatch or obj$nbatch.blen or obj$blen.nspac or obj$nspac.outfun or obj$outfun.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.
obj.
The initial state must satisfy obj(state, ...) > 0.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])Run the code above in your browser using DataLab