Description Usage Arguments Value See Also Examples
View source: R/metropolisHastings.R
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.
1 2 3 |
init |
An parameter vector with which to initialise the MCMC algorithm. |
logLik |
A function which takes a parameter (such as |
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 |
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 |
iters |
The number of MCMC iterations required (_after_ thinning). |
thin |
The required thinning factor. eg. only store every |
verb |
Boolean indicating whether some progress information
should be printed to the console. Defaults to |
debug |
Boolean indicating whether debugging information is required. Prints information about each iteration to console, to, eg., debug a crashing sampler. |
A matrix with rows representing samples from the posterior distribution.
pfMLLik
, StepGillespie
, abcRun
,
simTs
, stepLVc
, metrop
1 2 3 4 5 6 7 8 9 10 11 12 13 | ## 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)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.