Sample from a (log)density using hybrid Monte Carlo
Description
hybridMC()
samples from a (log) joint density function defined up to a proportionality constant. It uses the
multipoint hybrid Monte Carlo methods described by Liu (2001). The function supports a full range of tweaking options to
minimize autocorrlation.
Usage
1 
Arguments
y.start 
A vector of starting values 
n.samp 
The number of samples to draw. Each previous value is used as the starting value for the next sample 
logDens 
The logdensity function (up to an additive constant) from which you would like to sample. 
dLogDens 
The function giving the derivative of the logdensity function with respect to each variable. 
epsilon 
Either a single positive value giving the size of the time simulation steps, or a vector of length 2 giving the lower and upper bounds of the interval from which 
LFsteps 
An integer giving the number of leapfrog simulation steps to do 
compWeights 
The “masses” of the dimensions. Must be either a single numeric value or a vector of the same length as 
MPwidth 
The (integer) size of the multipoint window. The default is 1, meaning no multipoint

MPweights 
A vector of length 
progress 
An integer giving the number of samples between updates to a text progress bar. If 
... 
Arguments to pass to 
Details
The density should have support of the whole real line for every dimension. The quality of the samples can be improved by tweaking epsilon
, LFsteps
, MPwidth
, and to to a lesser extent, compWeights
and MPweights
.
If you use the progress bar, keep in mind that updating the progress bar takes a little time. Too many updates will slow down your sampling, so pick a reasonable value for progress
.
Value
hybridMC()
returns an object of type mcmc
(from the coda
package). Each row is a sample from the joint density.
Author(s)
Richard D. Morey
References
Liu (2001) "Monte Carlo strategies in scientific computing"
Examples
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22  #### Example 1: Jointly sample from two independent double exponentials
## log density function for double exponential
L = function(x)
sum(abs(x))
dL = function(x)
sign(x)
startVal=c(1,1)
samples = hybridMC(y.start=startVal, n.samp=5000, logDens=L, dLogDens=dL, epsilon=.2, LFsteps=10)
# Plot the MCMC chains and densities
plot(samples)
# Plot a histogram of the first variable, with true density
hist(samples[,1],freq=FALSE,breaks=50)
x = seq(5,5,len=100)
lines(x,.5*dexp(abs(x)),col="red") #true density in red
# Autocorrelation function
acf(samples)

Want to suggest features or report bugs for rdrr.io? Use the GitHub issue tracker. Vote for new features on Trello.