mh_sampler: Random Walk Metropolis-Hastings sampler using multiple...

In svdataman/tonic: Generate and test MCMC output using RW-MH or Goodman-Weare methods

Description

mh_sampler returns posterior samples from a RW-MH sampler.

Usage

mh_sampler(posterior, theta.0, nsteps = 10000, nchains = 5,
  burn.in = 2000, update = 5, chatter = 0, cov = NULL,
  thin = NULL, merge.chains = TRUE, adapt = FALSE,
  proposal = "normal", ...)

Arguments

posterior	(function) name of the log(densiy) function to sampling from
theta.0	(vector) initial values of the M variables
nsteps	(integer) total number of samples required
nchains	(integer) number of distinct 'chains' to run.
burn.in	(integer) the 'burn-in' period for the ensemble.
update	(integer) print a progress update after update steps
chatter	(integer) how verbose is the output? (0=quiet, 1=normal, 2=verbose)
cov	(array/matrix) covariance matrix for proposal distribution.
thin	(integer) keep only every <thin> sample
merge.chains	(logical) combine output from all chains into one
adapt	(logical) if TRUE then adjust the covariance matrix using the samples from the burn in period.
proposal	(string) Use "normal" proposal, or else "t" distribution.
...	(anything else) other arguments needed for posterior function

Details

A simple, pure R function to carry generate MCMC output using the random walk Metropolis-Hastings algorithm. Uses a Normal or t distribution for the proposal.

Value

A list with components

theta	(array) nsteps samples from M-dimensional posterior [nsteps rows, M columns]
func	(string) name of posterior function sampled
lpost	(vector) nsteps values of the LogPosterior density at each sample position
method	(string) sampling method used ("mh_sampler")
nchains	(integer) number of chains used
accept.rate	(float) the fraction of proposals accepted.
fcall	Full list of the function call.

If merge.chains = FALSE then theta will be a 3D array with dimensions nchains * (nsteps/nchains) * M.

Notes

Generate samples from some 'target' density function specified by the posterior function. Uses the random walk Metropolis-Hastings algorithm. The 'proposal' distribution is either a multivariate normal (default) or a Student's t distribution (with df=3). The latter has slightly fatter tails, so might be preferred for non-normal target densities.

This function generates a total of nsteps samples after a 'burn in' period of burn.in steps. Allowing for a burn in period helps remove memory of the starting position (important if this is not well chosen). We run nchains in parallel, from slightly different starting positions. Each chain only needs to run for nsteps / nchain iterations (after the burn in is complete). Having multiple chains helps assess convergence: are they 'well mixed' meanng that samples from each walker appear to be drawn from the same distribution? We can examine the moments and histograms of the chains to check this.

The proposal density (normal or t) requires a covariance matrix. This is specified on input by the cov argument. If not specified this will default to an identity matrix. Chosing a sensible covariance matrix for the proposal is part of the 'art' of MCMC. If the proposal is badly chosen - and nothing like the target density - then convergence (in distribution) will be very slow.

MCMC practitioners often aim for an acceptance rate of 0.1-0.5. If the acceptance rate is very low it means there will be many repeats of the same values before new ones are accepted, and so the random walk will explore the target density slowly. This is usually because the proposal density is too large compared to the target density - try reducing the scale of the covariance matrix. If the acceptance rate is very high (approaching 100%) then probably the proposal distribution is too small, each new proposed position is very close to the current position (and so of similar density and highly likely to be accepted). Again, this means the random walk will move only slowly through the target density. In this case, making the covariance larger usually helps.

The target density function: The target density should is specified by the posterior function. In fact, this function should compute the log density given parameters theta and any other arguments, i.e. log(p) = posterior(theta, ...) where the vector of parameters theta is the first argument of the posterior function. Where the density is zero or not defined, e.g. because prior = 0 for certain values of the parameters, it should return -Inf. Otherwise, the output of posterior(theta, ...) should be a real, scalar value.

See Also

chain_convergence, gw_sampler, chain_diagnosis, contour_matrix

Examples

my_posterior <- function(theta) {
  cov <- matrix(c(1,0.98,0.8,0.98,1.0,0.97,0.8,0.97,2.0), nrow = 3)
  logP <- mvtnorm::dmvnorm(theta, mean = c(-1, 2, 0), sigma = cov, log = TRUE)
  return(logP)
}
chain <- mh_sampler(my_posterior, theta.0 = c(0,0,0), nsteps = 10e4, burn.in = 1e4) 

svdataman/tonic documentation built on Aug. 2, 2019, 3:21 p.m.