R: The parallel Tempering algorithm

parallelTempering {EMC}

R Documentation

The parallel Tempering algorithm

Description

Given a multi-modal and multi-dimensional target density function, a (possibly asymmetric) proposal distribution and a temperature ladder, this function produces samples from the target using the parallel tempering algorithm.

Below sampDim refers to the dimension of the sample space, temperLadderLen refers to the length of the temperature ladder, and levelsSaveSampForLen refers to the length of the levelsSaveSampFor.

Usage

parallelTempering(nIters,                         
                  temperLadder,                   
                  startingVals,                   
                  logTarDensFunc,                 
                  MHPropNewFunc,                  
                  logMHPropDensFunc = NULL,       
                  MHBlocks          = NULL,
                  MHBlockNTimes     = NULL,
                  moveProbsList     = NULL,           
                  moveNTimesList    = NULL,          
                  levelsSaveSampFor = NULL,
                  nThin             = 1,
                  saveFitness       = FALSE,            
                  verboseLevel      = 0,
                  ...)

Arguments

`nIters`	`integer` > 0.
`temperLadder`	`double` vector with all positive entries, in decreasing order.
`startingVals`	`double` matrix of dimension `temperLadderLen` x `sampDim` or vector of length `sampDim`, in which case the same starting values are used for every temperature level.
`logTarDensFunc`	`function` of two arguments `(draw, ...)` that returns the target density evaluated in the log scale.
`MHPropNewFunc`	`function` of four arguments `(temperature, block, currentDraw, ...)` that returns new Metropolis-Hastings proposals. See details below on the argument block.
`logMHPropDensFunc`	`function` of five arguments `(temperature, block, currentDraw, proposalDraw, ...)` that returns the proposal density evaluated in the log scale. See details below on the argument block.
`MHBlocks`	`list` of integer vectors giving dimensions to be blocked together for sampling. It defaults to `as.list(1:sampDim)`, i.e., each dimension is treated as a block on its own. See details below for an example.
`MHBlockNTimes`	`integer` vector of number of times each block given by `MHBlocks` should be sampled in each iteration. It defaults to `rep(1, length(MHBlocks))`. See details below for an example.
`moveProbsList`	named `list` of probabilities adding upto 1.
`moveNTimesList`	named `list` of integers >= 0.
`levelsSaveSampFor`	`integer` vector with positive entries.
`nThin`	`integer` >= 1. Every `nThin` draw is saved.
`saveFitness`	`logical`.
`verboseLevel`	`integer`, a value >= 2 produces a lot of output.
`...`	optional arguments to be passed to `logTarDensFunc`, `MHPropNewFunc` and `logMHPropDensFunc`.

Details

MHPropNewFunc and logMHPropDensFunc: The MHPropNewFunc and the logMHPropDensFunc are called multiple times by varying the block argument over 1:length(MHBlocks), so these functions should know how to generate a proposal from the currentDraw or to evaluate the proposal density depending on which block was passed as the argument. See the example section for sample code.
MHBlocks and MHBlockNTimes: Blocking is an important and useful tool in MCMC that helps speed up sampling and hence mixing. Example: Let sampDim = 6. Let we want to sample dimensions 1, 2, 4 as one block, dimensions 3 and 5 as another and treat dimension 6 as the third block. Suppose we want to sample the three blocks mentioned above 1, 5 and 10 times in each iteration, respectively. Then we could set MHBlocks = list(c(1, 2, 4), c(3, 5), 6) and MHBlockNTimes = c(1, 5, 10)

The parallel tempering algorithm: The parallel tempering (PT; Liang and Wong, 2001) algorithm is composed of the following moves:

MH	Metropolis-Hastings or mutation
RE	(random) exchange

The current function could be used to run the PT algorithm by specifying what moves to employ using the following variables.

moveProbsList and moveNTimesList: The allowed names for components of moveProbsList and moveNTimesList come from the abbreviated names of the moves above. For example, the following specifications are valid:

moveProbsList = list(MH = 0.4,
                     RE = 0.6)

moveNTimesList = list(MH = 1,
                      RE = temperLadderLen)

levelsSaveSampFor

By default, samples are saved and returned for temperature level temperLadderLen. The levelsSaveSampFor could be used to save samples from other temperature levels as well (e.g., levelsSaveSampFor = 1:temperLadderLen saves samples from all levels).

saveFitness

The term fitness refers to the function H(x), where the target density of interest is given by:

g(x) \propto \exp[ -H(x) / τ_{min} ]

H(x) is also known as the energy function. By default, the fitness values are not saved, but one can do so by setting saveFitness = TRUE.

Value

Below nSave refers to ceil(nIters / nThin). This function returns a list with the following components:

`draws`	`array` of dimension `nSave` x `sampDim` x `levelsSaveSampForLen`, if `saveFitness = FALSE`. If `saveFitness = TRUE`, then the returned array is of dimension `nSave` x `(sampDim + 1)` x `levelsSaveSampForLen`; i.e., each of the `levelsSaveSampForLen` matrices contain the fitness values in their last column.
`acceptRatios`	`matrix` of the acceptance rates for various moves used.
`detailedAcceptRatios`	`list` of matrices with detailed summary of the acceptance rates for various moves used.
`nIters`	the `nIters` argument.
`nThin`	the `nThin` argument.
`nSave`	as defined above.
`temperLadder`	the `temperLadder` argument.
`startingVals`	the `startingVals` argument.
`moveProbsList`	the `moveProbsList` argument.
`moveNTimesList`	the `moveNTimesList` argument.
`levelsSaveSampFor`	the `levelsSaveSampFor` argument.
`time`	the time taken by the run.

Note

The effect of leaving the default value NULL for some of the arguments above are as follows:

`logMHPropDensFunc`	the proposal density `MHPropNewFunc` is deemed symmetric.
`MHBlocks`	`as.list(1:sampDim)`.
`MHBlockNTimes`	`rep(1, length(MHBlocks))`.
`moveProbsList`	`list(MH = 0.4, RC = 0.3, SC = 0.3)`.
`moveNTimesList`	`list(MH = 1, RC = mm, SC = mm, RE = nn)`, where
	`mm <- floor(nn / 2)` and `nn <- temperLadderLen`.
`levelsSaveSampFor`	`temperLadderLen`.

Author(s)

Gopi Goswami goswami@stat.harvard.edu

References

Faming Liang and Wing H.Wong (2001). Real-Parameter Evolutionary Monte Carlo with Applications to Bayesian Mixture Models. Journal of the American Statistical Association 96:653-666.

Examples

## Not run: 
samplerObj <-
    with(VShapedFuncGenerator(-13579),
         parallelTempering(nIters            = 2000,
                           temperLadder      = c(15, 6, 2, 1),
                           startingVals      = c(0, 0),
                           logTarDensFunc    = logTarDensFunc,
                           MHPropNewFunc     = MHPropNewFunc,
                           levelsSaveSampFor = seq_len(4), 
                           verboseLevel      = 1))
print(samplerObj)
print(names(samplerObj))
with(samplerObj,
 {
     print(detailedAcceptRatios)
     print(dim(draws))
     par(mfcol = c(2, 2))
     for (ii in seq_along(levelsSaveSampFor)) {
         main <- paste('temper:', round(temperLadder[levelsSaveSampFor[ii]], 3))
         plot(draws[ , , ii],
              xlim = c(-5, 20),
              ylim = c(-8, 8),
              pch  = '.',
              ask  = FALSE,
              main = as.expression(main),
              xlab = as.expression(substitute(x[xii], list(xii = 1))),
              ylab = as.expression(substitute(x[xii], list(xii = 2))))
     }
 })

## End(Not run)

[Package EMC version 1.3 Index]