EMC (version 1.0)

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

Usage

MetropolisHastings(nIters,
                   startingVal,         
                   logTarDensFunc,      
                   propNewFunc,         
                   logPropDensFunc,     
                   MHBlocks      = NULL,
                   MHBlockNTimes = NULL,
                   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.

saveFitness

logical indicating whether fitness values should be saved. See details below.

verboseLevel

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

...

optional arguments to be passed to logTarDensFunc, propNewFunc and logPropDensFunc.

Value

• This function returns a list with the following components:
• draws
  matrix of dimension nIters $\times$ sampDim, if saveFitness = FALSE. If saveFitness = TRUE, then the returned matrix is of dimension nIters $\times$ (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.
• startingVal
  the startingVal argument.
• time
  the time taken by the run.

Details

[object Object],[object Object],[object Object]

References

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

See Also

randomWalkMetropolis, parallelTempering, evolMonteCarlo

Examples

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

    means          <- c(1, 1)
    disp           <- ARDisp(-0.95)
    logTarDensFunc <-
        function (draw, ...)
            dmvnorm(draw, means, disp, log = TRUE)

    tDF         <- 3
    propNewFunc <-
        function (block, currentDraw, ...)
        {
            proposalDraw        <- currentDraw
            proposalDraw[block] <- rt(1, tDF)
            proposalDraw
        }            
    
    logPropDensFunc <-
        function (block, currentDraw, proposalDraw, ...)
            dt(proposalDraw[block], tDF, log = TRUE)

    list(logTarDensFunc  = logTarDensFunc,
         propNewFunc     = propNewFunc,
         logPropDensFunc = logPropDensFunc)
}


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