metropolisHastings {smfsb}

Run a Metropolis-Hastings MCMC algorithm for the parameters of a Bayesian posterior distribution

Description

Run a Metropolis-Hastings MCMC algorithm for the parameters of a Bayesian posterior distribution. Note that the algorithm carries over the old likelihood from the previous iteration, making it suitable for problems with expensive likelihoods, and also for "exact approximate" pseudo-marginal or particle marginal MH algorithms.

Usage

metropolisHastings(init, logLik, rprop, dprop=function(new, old, ...){1},
                   dprior=function(x, ...){1}, iters=10000, thin=10,
                   verb=TRUE, debug=FALSE)

Arguments

init

An parameter vector with which to initialise the MCMC algorithm.

logLik

A function which takes a parameter (such as init) as its only required argument and returns the log-likelihood of the data. Note that it is fine for this to return the log of an unbiased estimate of the likelihood, in which case the algorithm will be an "exact approximate" pseudo-marginal MH algorithm.

rprop

A function which takes a parameter as its only required argument and returns a single sample from a proposal distribution.

dprop

A function which takes a new and old parameter as its first two required arguments and returns the (log) density of the new value conditional on the old. It should accept an optional parameter log, and must return the log-density when log is TRUE. Defaults to a flat function which causes this term to drop out of the acceptance probability. It is fine to use the default for _any_ _symmetric_ proposal, since the term will also drop out for any symmetric proposal.

dprior

A function which take a parameter as its only required argument and returns the (log) density of the parameter value under the prior. It should accept an optional parameter log, and must return the log-density when log is TRUE. Defaults to a flat function which causes this term to drop out of the acceptance probability. People often use a flat prior when they are trying to be "uninformative" or "objective", but this is slightly naive. In particular, what is "flat" is clearly dependent on the parametrisation of the model.

iters

The number of MCMC iterations required (_after_ thinning).

thin

The required thinning factor. eg. only store every thin iterations.

verb

Boolean indicating whether some progress information should be printed to the console. Defaults to TRUE.

debug

Boolean indicating whether debugging information is required. Prints information about each iteration to console, to, eg., debug a crashing sampler.

Value

A matrix with rows representing samples from the posterior distribution.

See Also

pfMLLik, StepGillespie, abcRun, simTs, stepLVc, metrop

Examples

## First simulate some synthetic data
data = rnorm(250,5,2)
## Now use MH to recover the parameters
llik = function(x) { sum(dnorm(data,x[1],x[2],log=TRUE)) }
prop = function(x) { rnorm(2,x,0.1) }
prior = function(x, log=TRUE) {
    l = dnorm(x[1],0,100,log=TRUE) + dgamma(x[2],1,0.0001,log=TRUE)
    if (log) l else exp(l)
}
out = metropolisHastings(c(mu=1,sig=1), llik, prop,
                         dprior=prior, verb=FALSE)
out = out[1000:10000,]
mcmcSummary(out, truth=c(5,2), rows=2, plot=FALSE)