harmonicHMC: Sample from a truncated Gaussian distribution with the harmonic HMC

Description

Generate MCMC samples from a d-dimensional truncated Gaussian distribution with constraints Fx+g >= 0 using the Harmonic Hamiltonian Monte Carlo sampler (Harmonic-HMC).

Usage

harmonicHMC(
  nSample,
  burnin = 0,
  mean,
  choleskyFactor,
  constrainDirec,
  constrainBound,
  init,
  time = c(pi/8, pi/2),
  precFlg,
  seed = 1,
  extraOutputs = c()
)

Value

When extraOutputs is empty (default), returns an nSample-by-d matrix of samples.

When extraOutputs contains "numBounces" and/or "bounceDistances", returns a list with elements:

samples: nSample-by-d matrix of samples
numBounces: Vector of bounce counts per sample (if requested)
bounceDistances: List of bounce distances per sample (if requested)

Arguments

nSample: number of samples after burn-in.
burnin: number of burn-in samples (default = 0).
mean: a d-dimensional mean vector.
choleskyFactor: upper triangular matrix R from Cholesky decomposition of precision or covariance matrix into R^TR.
constrainDirec: the k-by-d F matrix (k is the number of linear constraints).
constrainBound: the k-dimensional g vector.
init: a d-dimensional vector of the initial value. init must satisfy all constraints.
time: HMC integration time for each iteration. Can either be a scalar value for a fixed time across all samples, or a length 2 vector of a lower and upper bound for uniform distribution from which the time is drawn from for each iteration.
precFlg: logical. whether choleskyFactor is from precision (TRUE) or covariance matrix (FALSE).
seed: random seed (default = 1).
extraOutputs: vector of strings. "numBounces" and/or "bounceDistances" can be requested, with the latter containing the distances in-between bounces for each sample and hence incurring significant computational and memory costs.

References

Pakman, A. and Paninski, L. (2014). Exact Hamiltonian Monte Carlo for Truncated Multivariate Gaussians. Journal of Computational and Graphical Statistics. doi:10.1080/10618600.2013.788448

Examples

Run this code

set.seed(1)
d <- 10
A <- matrix(runif(d^2)*2 - 1, ncol=d)
precMat <- t(A) %*% A
R <- cholesky(precMat)
mu <- rep(0, d)
constrainDirec <- diag(d)
constrainBound <- rep(0,d)
initial <- rep(1, d)
results <- harmonicHMC(1000, 1000, mu, R, constrainDirec, constrainBound, initial, precFlg = TRUE)