EMC (version 1.0)

parallelTempering: 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.

Usage

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

Arguments

nIters

integer $>$ 0.

temperLadder

double vector with all positive entries, in decreasing order.

startingVals

double matrix of dimension temperLadderLen $\times$ 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 $\ge$ 0.

levelsSaveSampFor

integer vector with positive entries.

saveFitness

logical.

verboseLevel

integer, a value $\ge$ 2 produces a lot of output.

...

optional arguments to be passed to logTarDensFunc, MHPropNewFunc and logMHPropDensFunc.

Value

• This function returns a list with the following components:
• draws
  array of dimension nIters $\times$ sampDim $\times$ levelsSaveSampForLen, if saveFitness = FALSE. If saveFitness = TRUE, then the returned array is of dimension nIters $\times$ (sampDim + 1) $\times$ 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.
• 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.

Details

[object Object],[object Object],.,[object Object],rl{ 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.

[object Object],moveProbsList = list(MH = 0.4, RE = 0.6)

moveNTimesList = list(MH = 1, RE = temperLadderLen),[object Object],[object Object]

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.

See Also

evolMonteCarlo

Examples

## The V-shaped distribution
VShapedFuncGenerator <-
    function (seed = 13579)
{
    set.seed(seed)    
    dd     <- 2
    ARDisp <-
        function (rho)
        {
            tmp <- rep(1, dd)
            diag((1 - rho) * tmp) + rho * tmp %*% t(tmp)
        }

    nMixComps  <- 2
    logWeights <- log(rep(1 / nMixComps, nMixComps))
    meanMat    <- matrix(c(1, 1, 15, 1), nMixComps, byrow = TRUE)
    dispParam  <- c(-0.95, 0.95)
    dispArr    <- array(dim = c(2, 2, nMixComps))
    for (ii in seq_len(nMixComps)) {
        dispArr[ , , ii] <- ARDisp(dispParam[ii])
    }
    logTarDensFunc <-
        function (draw, ...)
        {
            ld <- sapply(seq_len(nMixComps), FUN =
                         function (ii)
                     {
                         dmvnorm(draw, meanMat[ii, ], dispArr[ , , ii], log = TRUE)
                     })
            ww <- logWeights + ld
            mm <- max(ww)
            mm + log(sum(exp(ww - mm)))
        }
    
    MHProposalSD  <- c(1.0, 1.0)
    MHPropNewFunc <-
        function (temperature, block, currentDraw, ...)
        {
            proposalDraw        <- currentDraw
            proposalDraw[block] <- rnorm(1, currentDraw[block],
                                         sqrt(temperature) * MHProposalSD[block])
            proposalDraw
        }

    list(logTarDensFunc = logTarDensFunc,
         MHPropNewFunc  = MHPropNewFunc)
}


samplerObj <-
    with(VShapedFuncGenerator( ),
         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))))
     }
 })