Simulated Tempering and Umbrella Sampling

Description

Markov chain Monte Carlo for continuous random vectors using parallel or serial simulated tempering, also called umbrella sampling. For serial tempering the state of the Markov chain is a pair (i, x), where i is an integer between 1 and k and x is a vector of length p. This pair is represented as a single real vector c(i, x). For parallel tempering the state of the Markov chain is vector of vectors (x[1], …, x[k]), where each x is of length p. This vector of vectors is represented as a k by p matrix.

Usage

temper(obj, initial, neighbors, nbatch, blen = 1, nspac = 1, scale = 1,
    outfun, debug = FALSE, parallel = FALSE, ...)

Arguments

Details

Serial tempering simulates a mixture of distributions of a continuous random vector. The number of components of the mixture is 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 to k] h(i, x)

Parallel tempering simulates a product of distributions of a continuous random vector. Denote the state (x[1], …, x[k]), then the unnormalized joint density of the equilibrium distribution is

h(x[1], …, x[k]) = prod[i = 1 to 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 “random-walk” Metropolis updates with multivariate normal proposal, the proposal distribution being determined by the argument 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 preserving h(x[1], …, 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.

Value

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

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

Warning

If 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.

Examples

d <- 9
witch.which <- c(0.1, 0.3, 0.5, 0.7, 1.0)
ncomp <- length(witch.which)

neighbors <- matrix(FALSE, ncomp, ncomp)
neighbors[row(neighbors) == col(neighbors) + 1] <- TRUE
neighbors[row(neighbors) == col(neighbors) - 1] <- TRUE

ludfun <- function(state, log.pseudo.prior = rep(0, ncomp)) {
    stopifnot(is.numeric(state))
    stopifnot(length(state) == d + 1)
    icomp <- state[1]
    stopifnot(icomp == as.integer(icomp))
    stopifnot(1 <= icomp && icomp <= ncomp)
    stopifnot(is.numeric(log.pseudo.prior))
    stopifnot(length(log.pseudo.prior) == ncomp)
    theta <- state[-1]
    if (any(theta > 1.0)) return(-Inf)
    bnd <- witch.which[icomp]
    lpp <- log.pseudo.prior[icomp]
    if (any(theta > bnd)) return(lpp)
    return(- d * log(bnd) + lpp)
}

# parallel tempering
thetas <- matrix(0.5, ncomp, d)
out <- temper(ludfun, initial = thetas, neighbors = neighbors, nbatch = 20,
    blen = 10, nspac = 5, scale = 0.56789, parallel = TRUE, debug = TRUE)

# serial tempering
theta.initial <- c(1, rep(0.5, d))
# log pseudo prior found by trial and error
qux <- c(0, 9.179, 13.73, 16.71, 20.56)

out <- temper(ludfun, initial = theta.initial, neighbors = neighbors,
    nbatch = 50, blen = 30, nspac = 2, scale = 0.56789,
    parallel = FALSE, debug = FALSE, log.pseudo.prior = qux)