c(i, x). For parallel tempering the state of the Markov chain
is vector of vectors $(x_1, \ldots, x_k)$,
where each x is of length $p$. This vector of vectors is
represented as a $k \times p$ matrix.temper(obj, initial, neighbors, nbatch, blen = 1, nspac = 1, scale = 1,
outfun, debug = FALSE, parallel = FALSE, ...)"tempering" from
a previous run. If a function, should evalutate the log unnormalized
density $\log h(i, x)$ of the desired equilibrium
distribution of the Markov chain for serial temperic(i, x) as
described above. For parallel tempering, a real
$k \times p$ matrix as described above. In either case,
the initial state of the Markov chain.k
by k. Elements that are TRUE indicate jumps
or swaps attempted by the Markov chain.outfun(state, ...) are returned. The argument state
is like the argument initial of this function. If missing, the
batch means of the rTRUE extra output useful for testing.TRUE does parallel tempering, if FALSE does
serial tempering.obj or outfun."mcmc", subclass "tempering",
which is a list containing at least the following components:outfun is missing, an nbatch by k by p
array. Otherwise, an nbatch by m matrix, where m
is the length of the result of outfun.nbatch
by k matrix giving batch means for the multivariate Bernoulli
random vector that is all zeros except for a 1 in the i-th place
when the current state is $(i, x)$.k giving the acceptance rate for each component.k by k matrix giving the acceptance rate for each allowed
jump or swap component. NA for elements such that the corresponding
elements of neighbors is FALSE.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, then the log unnormalized
density function 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 temper and that works too.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 temper,
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.
k, and the dimension
of the random vector is p. Denote the state $(i, x)$, where $i$
is a positive integer between 1 and $k$, and let $h(i, x)$ denote
the unnormalized joint density of their equilibrium distribution.
The logarithm of this function is what obj or obj$lud calculates.
The mixture distribution is the marginal for $x$ derived from
the equilibrium distribution $h(i, x)$, that is,
$$h(x) = \sum_{i = 1}^k h(i, x)$$Parallel tempering simulates a product of distributions of a continuous random vector. Denote the state $(x_1, \ldots, x_k)$, then the unnormalized joint density of the equilibrium distribution is $$h(x_1, \ldots, x_k) = \prod_{i = 1}^k h(i, x_i)$$
The update mechanism of the Markov chain combines two kinds of elementary updates: jump/swap updates (jump for serial tempering, swap for parallel tempering) and within-component updates. Each iteration of the Markov chain one of these elementary updates is done. With probability 1/2 a jump/swap update is done, and with probability 1/2 a with-component update is done.
Within-component updates are the same for both serial and parallel tempering.
They are scale. In serial tempering, the $x$ part of the current state
$(i, x)$ is updated preserving $h(i, x)$.
In parallel tempering, an index $i$ is chosen at random and the part
of the state $x_i$ representing that component is updated,
again preserving $h(i, x)$.
Jump updates choose uniformly at random a neighbor of the current component:
if $i$ indexes the current component, then it chooses uniformly at random
a $j$ such that neighbors[i, j] == TRUE. It then does does a
Metropolis-Hastings update for changing the current state from $(i, x)$
to $(j, x)$.
Swap updates choose a component uniformly at random and a neighbor of that
component uniformly at random: first an index $i$ is chosen uniformly
at random between 1 and $k$, then an index $j$ is chosen
uniformly at random such that neighbors[i, j] == TRUE. It then does
does a Metropolis-Hastings update for swapping the states of the
two components: interchanging $x_i$ and $x_j$
while perserving $h(x_1, \ldots, x_k)$.
The initial state must satisfy lud(initial, ...) > - Inf for serial
tempering or must satisfy lud(initial[i, ], ...) > - Inf for each
i for parallel tempering, where lud is either obj
or obj$lud.
That is, the initial state must have positive probability.
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