evolMonteCarlo: evolutionary Monte Carlo algorithm

In EMC: Evolutionary Monte Carlo (EMC) 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 evolutionary Monte Carlo 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

evolMonteCarlo(nIters,
               temperLadder,                       
               startingVals,                       
               logTarDensFunc,                     
               MHPropNewFunc,                      
               logMHPropDensFunc = NULL,           
               MHBlocks          = NULL,
               MHBlockNTimes     = NULL,
               moveProbsList     = NULL,               
               moveNTimesList    = NULL,              
               SCRWMNTimes       = NULL,                     
               SCRWMPropSD       = 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.
SCRWMNTimes	integer > 0.
SCRWMPropSD	double > 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 EMC and the TOEMC algorithm

The evolutionary Monte Carlo (EMC; Liang and Wong, 2001) algorithm is composed of the following moves:

MH	Metropolis-Hastings or mutation
RC	real crossover
SC	snooker crossover
RE	(random) exchange

The target oriented EMC (TOEMC; Goswami and Liu, 2007) algorithm has the following additional moves on top of EMC:

BCE	best chromosome exchange
BIRE	best importance ratio exchange
BSE	best swap exchange
CE	cyclic exchange

The current function could be used to run both EMC and TOEMC algorithms 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,
                     RC = 0.3,
                     SC = 0.3)
      

moveNTimesList = list(MH = 1, 
                      RC = floor(temperLadderLen / 2),
                      SC = floor(temperLadderLen / 2),
                      RE = temperLadderLen)           
      

SCRWMNTimes and SCRWMPropSD

The conditional sampling step of the snooker crossover (SC) move is done using random walk Metropolis (RWM) with normal proposals; SCRWMNTimes and SCRWMPropSD are the number of RWM draws and the proposal standard deviation for the RWM step, respectively. Note these variables are only required if the SC move has positive probability in moveProbsList or a positive number of times in moveNTimesList.

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.
SCRWMNTimes	1, if SC is used.
SCRWMPropSD	needs to be provided by the user, if SC is used.
levelsSaveSampFor	temperLadderLen.

Author(s)

Gopi Goswami [e-mail]

References

Gopi Goswami and Jun S. Liu (2007). On learning strategies for evolutionary Monte Carlo. Statistics and Computing 17:1:23-38.

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

parallelTempering

Examples

## Not run: 
samplerObj <-
    with(VShapedFuncGenerator(-13579),
     {
         allMovesNTimesList <- list(MH = 1, RC = 2, SC = 2, RE = 4,
                                    BIRE = 2, BCE = 2, BSE = 2, CE = 2)
         evolMonteCarlo(nIters            = 2000,
                        temperLadder      = c(15, 6, 2, 1),
                        startingVals      = c(0, 0),
                        logTarDensFunc    = logTarDensFunc,
                        MHPropNewFunc     = MHPropNewFunc,
                        moveNTimesList    = allMovesNTimesList,
                        SCRWMNTimes       = 1,
                        SCRWMPropSD       = 3.0,
                        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)

Example output

Loading required package: mvtnorm
Loading required package: MASS
##
## Evolutionary Monte Carlo Package (EMC)
##
## Functionality: random walk Metropolis, general Metropolis-Hastings
## parallel tempering, target oriented EMC (TOEMC), temperature ladder
## construction and placement
##
## Use: "help(package = EMC)" at the R prompt for more info
##
## Copyright (C) 2006-2018 Gopi Goswami
##
##    Created by: Gopi Goswami <[e-mail]>
## Maintained by: Gopi Goswami <[e-mail]>
##


BEGIN: EMC
..........
[Time to finish (est):        3 secs, this iter:        100]....................
[Time to finish (est):        3 secs, this iter:        300]....................
[Time to finish (est):        2 secs, this iter:        500]....................
[Time to finish (est):        2 secs, this iter:        700]....................
[Time to finish (est):        2 secs, this iter:        900]....................
[Time to finish (est):        2 secs, this iter:       1100]....................
[Time to finish (est):        1 secs, this iter:       1300]....................
[Time to finish (est):        1 secs, this iter:       1500]....................
[Time to finish (est):        1 secs, this iter:       1700]....................
[Time to finish (est):        1 secs, this iter:       1900]..........
[Total time:          2 secs (usr),     0 secs (sys), this iter:       2000]
E N D: EMC

The temperature ladder:
[1] 15  6  2  1

The overall acceptance rate summary:
         ratio accepted proposed moveNTimes
MH   0.3754759     2367     6304          1
RC   0.2335526      284     1216          2
SC   0.3865894      467     1208          2
RE   0.5192500     4154     8000          4
BCE  0.8492500     3397     4000          2
BIRE 0.9437500     3775     4000          2
BSE  0.4750000     1900     4000          2
CE   0.2762500     1105     4000          2

 [1] "draws"                "acceptRatios"         "nIters"              
 [4] "nThin"                "nSave"                "levelsSaveSampFor"   
 [7] "temperLadder"         "startingVals"         "moveProbsList"       
[10] "moveNTimesList"       "detailedAcceptRatios" "time"                
$MH
       argminBlock   min argmaxBlock   max
level1           2 0.410           1 0.444
level2           2 0.365           1 0.367
level3           1 0.341           2 0.376
level4           1 0.327           2 0.373

$BetweenAllLevels
     argminLevels   min argmaxLevels   max
RC      4 <-->  3 0.159    1 <-->  2 1.000
SC      4 <-->  2 0.316    1 <-->  2 0.593
RE      3 <-->  2 0.396    3 <-->  4 0.567
BCE     4 <-->  1 0.793    3 <-->  4 0.879
BIRE    3 <-->  2 0.919    1 <-->  4 0.981
BSE     4 <-->  1 0.398    3 <-->  4 0.565

$TargetWithOtherLevels
     argminLevels   min argmaxLevels   max
RC      4 <-->  3 0.159    4 <-->  1 0.386
SC      4 <-->  2 0.316    4 <-->  1 0.518
RE      4 <-->  3 0.567    4 <-->  3 0.567
BCE     4 <-->  1 0.793    4 <-->  3 0.879
BIRE    4 <-->  2 0.958    4 <-->  1 0.981
BSE     4 <-->  1 0.398    4 <-->  3 0.565

$CE
              argminSelLength   min argmaxSelLength   max
ladderLength3               2 0.264               1 0.532
ladderLength4               3 0.124               2 0.493

[1] 2000    2    4

EMC documentation built on May 2, 2019, 5:11 a.m.