R: The Metropolis-Hastings algorithm

MetropolisHastings {EMC}

R Documentation

The Metropolis-Hastings algorithm

Description

Given a target density function and an asymmetric proposal distribution, this function produces samples from the target using the Metropolis Hastings algorithm.

Below sampDim refers to the dimension of the sample space.

Usage

MetropolisHastings(nIters,
                   startingVal,         
                   logTarDensFunc,      
                   propNewFunc,         
                   logPropDensFunc,     
                   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.
`logPropDensFunc`	`function` of four arguments `(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.
`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`, `propNewFunc` and `logPropDensFunc`.

Details

propNewFunc and logPropDensFunc: The propNewFunc and the logPropDensFunc 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)
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 goswami@stat.harvard.edu

References

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

Examples

## Not run: 
samplerObj <-
    with(CigarShapedFuncGenerator2(-13579),
         MetropolisHastings(nIters          = 5000,
                            startingVal     = c(0, 0),
                            logTarDensFunc  = logTarDensFunc,
                            propNewFunc     = propNewFunc,
                            logPropDensFunc = logPropDensFunc,
                            verboseLevel    = 2))
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))))    
 })

## End(Not run)

[Package EMC version 1.3 Index]