hmc_sampler: Hamiltonian MC sampler

In andrewzm/hmc: Hamiltonian Markov chain Monte Carlo

Description

Hamiltonian MC sampler

Usage

hmc_sampler(U, dUdq, M, eps_gen, L)

Arguments

U	potential energy function.
dUdq	derivative of potential energy function with respect to position.
M	the mass matrix (i.e.,~the covariance matrix in the kinetic energy function)
eps_gen	a function which generates the leapfrog step-size ε
L	number of leapfrog steps per proposal

Examples

### Model setup -- 2d correlated Gaussian
S <- matrix(c(1,-0.98,-0.98,1),2,2)
n <- nrow(S)
Q <- chol2inv(chol(S))
cholQ <- chol(Q)
U <- function(q) 0.5 * crossprod(cholQ %*% q)
dUdq <- function(q) Q %*% q
M <- diag(1,n)

### Sampler parameters -- set eps and L according to eigenvalues of covariance matrix
E <- eigen(S)
eps_gen <- function() round(min(E$values),2)
L = as.integer(max(E$values)/eps_gen())
print(paste0("eps = ",eps_gen(),". L = ",L))
sampler <- hmc_sampler(U = U,dUdq = dUdq, M = M,eps_gen = eps_gen,L = L)

### Now sample
N <- 1000
q <- matrix(0,n,N)
for(i in 2:N) q[,i] <- sampler(q = q[,(i-1)])
plot(t(q),ylim = c(-4,4),xlim=c(-4,4))

andrewzm/hmc documentation built on May 10, 2019, 11:15 a.m.