randomWalkMetropolis: The Random Walk Metropolis algorithm

In EMC: Evolutionary Monte Carlo (EMC) algorithm

Description

Given a target density function and a symmetric proposal generating function, this function produces samples from the target using the random walk Metropolis algorithm.

Below sampDim refers to the dimension of the sample space.

Usage

randomWalkMetropolis(nIters,              
                     startingVal,         
                     logTarDensFunc,      
                     propNewFunc,         
                     MHBlocks      = NULL,
                     MHBlockNTimes = NULL,
                     nThin         = 1,
                     saveFitness   = FALSE,
                     verboseLevel  = 0,
                     ...)    

Arguments

nIters	integer > 0.
startingVal	double vector of length sampDim.
logTarDensFunc	function of two arguments (draw, ...) that returns the target density evaluated in the log scale.
propNewFunc	function of three arguments (block, currentDraw, ...) that returns new Metropolis-Hastings proposals. 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.
nThin	integer >= 1. Every nThin draw is saved.
saveFitness	logical indicating whether fitness values should be saved. See details below.
verboseLevel	integer, a value >= 2 produces a lot of output.
...	optional arguments to be passed to logTarDensFunc and propNewFunc.

Details

propNewFunc

The propNewFunc is called multiple times by varying the block argument over 1:length(MHBlocks), so this function should know how to generate a proposal from the currentDraw 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)

saveFitness

The term fitness refers to the negative of the logTarDensFunc values. 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	matrix of dimension nSave x sampDim, if saveFitness = FALSE. If saveFitness = TRUE, then the returned matrix is of dimension nSave x (sampDim + 1), where the fitness values appear in its last column.
acceptRatios	matrix of the acceptance rates.
detailedAcceptRatios	matrix with detailed summary of the acceptance rates.
nIters	the nIters argument.
nThin	the nThin argument.
nSave	as defined above.
startingVal	the startingVal 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:

MHBlocks	as.list(1:sampDim).
MHBlockNTimes	rep(1, length(MHBlocks)).

Author(s)

Gopi Goswami [e-mail]

References

Jun S. Liu (2001). Monte Carlo strategies for scientific computing. Springer.

See Also

MetropolisHastings, parallelTempering, evolMonteCarlo

Examples

## Not run: 
samplerObj <-
    with(CigarShapedFuncGenerator1(-13579),
         randomWalkMetropolis(nIters         = 5000,
                              startingVal    = c(0, 0),
                              logTarDensFunc = logTarDensFunc,
                              propNewFunc    = propNewFunc,
                              verboseLevel   = 1))
print(samplerObj)
print(names(samplerObj))
with(samplerObj,
 {
     print(detailedAcceptRatios)
     print(dim(draws))
     plot(draws,
          xlim = c(-3, 5),
          ylim = c(-3, 4),
          pch  = '.',
          ask  = FALSE,
          main = as.expression(paste('# draws:', nIters)),
          xlab = as.expression(substitute(x[xii], list(xii = 1))),
          ylab = as.expression(substitute(x[xii], list(xii = 2))))    
 })


samplerObj <-
    with(threeDimNormalFuncGenerator(-13579),
     {
         randomWalkMetropolis(nIters          = 5000,
                              startingVal     = c(0, 0, 0),         
                              logTarDensFunc  = logTarDensFunc,     
                              propNewFunc     = propNewFunc,        
                              MHBlocks        = list(c(1, 2), 3),   
                              verboseLevel    = 1)                  
     })
print(samplerObj)
print(names(samplerObj))
with(samplerObj,
 {
     print(detailedAcceptRatios)
     print(dim(draws))
     pairs(draws,
           pch    = '.',                                                    
           ask    = FALSE,                                                  
           main = as.expression(paste('# draws:', nIters)),
           labels = c(as.expression(substitute(x[xii], list(xii = 1))), 
                      as.expression(substitute(x[xii], list(xii = 2))), 
                      as.expression(substitute(x[xii], list(xii = 3)))))
 })

## 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-2019 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):        2 secs, this iter:        600]....................
[Time to finish (est):        1 secs, this iter:       1100]....................
[Time to finish (est):        1 secs, this iter:       1600]....................
[Time to finish (est):        1 secs, this iter:       2100]....................
[Time to finish (est):        1 secs, this iter:       2600]....................
[Time to finish (est):        1 secs, this iter:       3100]....................
[Time to finish (est):        1 secs, this iter:       3600]....................
[Time to finish (est):        1 secs, this iter:       4100]....................
[Time to finish (est):        1 secs, this iter:       4600]................
[Total time:          1 secs (usr),     0 secs (sys), this iter:       5000]
E N D: EMC

The overall acceptance rate summary:
         ratio accepted proposed
overall 0.2719     2719    10000

[1] "draws"                "acceptRatios"         "nIters"              
[4] "nThin"                "nSave"                "detailedAcceptRatios"
[7] "time"                 "startingVal"         
             argminBlock   min argmaxBlock   max
acceptRatios           2 0.189           1 0.355
[1] 5000    2

BEGIN: EMC
..........
[Time to finish (est):        2 secs, this iter:        100]....................
[Time to finish (est):        2 secs, this iter:        600]....................
[Time to finish (est):        2 secs, this iter:       1100]....................
[Time to finish (est):        2 secs, this iter:       1600]....................
[Time to finish (est):        1 secs, this iter:       2100]....................
[Time to finish (est):        1 secs, this iter:       2600]....................
[Time to finish (est):        1 secs, this iter:       3100]....................
[Time to finish (est):        1 secs, this iter:       3600]....................
[Time to finish (est):        1 secs, this iter:       4100]....................
[Time to finish (est):        1 secs, this iter:       4600]................
[Total time:          2 secs (usr),     0 secs (sys), this iter:       5000]
E N D: EMC

The overall acceptance rate summary:
         ratio accepted proposed
overall 0.3279     3279    10000

[1] "draws"                "acceptRatios"         "nIters"              
[4] "nThin"                "nSave"                "detailedAcceptRatios"
[7] "time"                 "startingVal"         
             argminBlock   min argmaxBlock   max
acceptRatios           1 0.167           2 0.489
[1] 5000    3

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