harmonicHMC: Sample from a truncated Gaussian distribution with the...

In hdtg: Generate Samples from Multivariate Truncated Normal Distributions

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()
)

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.

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)

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

See Also

getHarmonicSample(), cholesky(), getInitialPosition()

Examples

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)

hdtg documentation built on Feb. 11, 2026, 5:07 p.m.