Efficient Cholesky decomposition

Description

Compute Cholesky decomposition of a matrix.

Usage

cholesky(A)

Arguments

Value

upper triangular matrix R such that A = U'U.

Create a Zigzag-HMC engine object

Description

Create the C++ object to set up SIMD vectorization for speeding up calculations for Zigzag-HMC ("Zigzag-HMC engine").

Usage

createEngine(
  dimension,
  lowerBounds,
  upperBounds,
  seed,
  mean,
  precision,
  flags = 128L
)

Arguments

Value

a list whose only element is the Zigzag-HMC engine object.

Create a Zigzag-NUTS engine object

Description

Create the C++ object to set up SIMD vectorization for speeding up calculations for Zigzag-NUTS ("Zigzag-NUTS engine").

Usage

createNutsEngine(
  dimension,
  lowerBounds,
  upperBounds,
  seed,
  stepSize,
  mean,
  precision,
  flags = 128L
)

Arguments

Value

a list whose only element is the Zigzag-NUTS engine object.

Draw a random Laplace momentum

Description

Generate a d-dimensional momentum where the density of each element is proportional to exp(-|pi|).

Usage

drawLaplaceMomentum(d)

Arguments

Value

a d-dimensional Laplace-distributed momentum.

Get an eligible initial value for a MTN with given mean and truncations

Description

For a given MTN the function returns an initial vector whose elements are one of: (1) middle point of the truncation interval if both lower and upper bounds are finite (2) lower (upper) bound +0.1 (-0.1) if only the lower (upper) bound is finite (3) the corresponding mean value if lower bound = -Inf are upper bound = Inf.

Usage

getInitialPosition(mean, lowerBounds, upperBounds)

Arguments

Value

an eligible d-dimensional initial vector.

Draw one Markovian zigzag sample

Description

Simulate the Markovian zigzag dynamics for a given position over a specified travel time.

Usage

getMarkovianZigzagSample(position, velocity = NULL, engine, travelTime)

Arguments

Value

A list containing the position and velocity after simulating the dynamics.

Draw one MTN sample with Zigzag-HMC or Zigzag-NUTS

Description

Simulate the Zigzag-HMC or Zigzag-NUTS dynamics on a given MTN.

Usage

getZigzagSample(position, momentum = NULL, nutsFlg, engine, stepZZHMC = NULL)

Arguments

Value

one MCMC sample from the target MTN.

Note

getZigzagSample is particularly efficient when the target MTN has a random mean and covariance/precision where one can reuse the Zigzag-HMC engine object while updating the mean and covariance. The following example demonstrates such a use.

Examples

set.seed(1)
n <- 1000
d <- 10
samples <- array(0, c(n, d))

# initialize MTN mean and precision
m <- rnorm(d, 0, 1)
prec <- rWishart(n = 1, df = d, Sigma = diag(d))[, , 1]
# call createEngine once
engine <- createEngine(dimension = d, lowerBounds = rep(0, d),
 upperBounds = rep(Inf, d), seed = 1, mean = m, precision = prec)

HZZtime <- sqrt(2) / sqrt(min(mgcv::slanczos(
 A = prec, k = 1,
 kl = 1
)[['values']]))

currentSample <- rep(0.1, d)
for (i in 1:n) {
  m <- rnorm(d, 0, 1)
  prec <- rWishart(n = 1, df = d, Sigma = diag(d))[,,1]
  setMean(sexp = engine$engine, mean = m)
  setPrecision(sexp = engine$engine, precision = prec)
  currentSample <- getZigzagSample(position = currentSample, nutsFlg = FALSE,
      engine = engine, stepZZHMC = HZZtime)
  samples[i,] <- currentSample
}

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 = NULL,
  extraOutputs = c()
)

Arguments

Value

samples: nSample-by-d matrix of samples or, if extraOutputs is non-empty, a list of samples and the extra outputs.

References

Pakman A, Paninski L (2014). “Exact Hamiltonian Monte Carlo for truncated multivariate Gaussians.” Journal of Computational and Graphical Statistics, 23(2), 518–542.

Examples

set.seed(1)
d <- 10
A <- matrix(runif(d^2)*2 - 1, ncol=d)
Sigma <- t(A) %*% A
R <- cholesky(Sigma)
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 = FALSE)

Set the mean for the target MTN

Description

Set the mean vector for a given Zigzag-HMC engine object.

Usage

setMean(sexp, mean)

Arguments

Set the precision matrix for the target MTN

Description

Set the precision matrix for a given Zigzag-HMC engine object.

Usage

setPrecision(sexp, precision)

Arguments

Sample from a truncated Gaussian distribution

Description

Generate MCMC samples from a d-dimensional truncated Gaussian distribution with element-wise truncations using the Zigzag Hamiltonian Monte Carlo sampler (Zigzag-HMC).

Usage

zigzagHMC(
  nSample,
  burnin = 0,
  mean,
  prec,
  lowerBounds,
  upperBounds,
  init = NULL,
  stepsize = NULL,
  nutsFlg = FALSE,
  precondition = FALSE,
  seed = NULL,
  diagnosticMode = FALSE
)

Arguments

Value

an nSample-by-d matrix of samples. If diagnosticMode is TRUE, a list with additional diagnostic information is returned.

References

Nishimura A, Zhang Z, Suchard MA (2024). “Zigzag path connects two Monte Carlo samplers: Hamiltonian counterpart to a piecewise deterministic Markov process.” Journal of the American Statistical Association, 1–13.

Nishimura A, Dunson DB, Lu J (2020). “Discontinuous Hamiltonian Monte Carlo for discrete parameters and discontinuous likelihoods.” Biometrika, 107(2), 365–380.

Examples

set.seed(1)
d <- 10
A <- matrix(runif(d^2)*2-1, ncol=d)
covMat <- t(A) %*% A
precMat <- solve(covMat)
initial <- rep(1, d)
results <- zigzagHMC(nSample = 1000, burnin = 1000, mean = rep(0, d), prec = precMat,
lowerBounds = rep(0, d), upperBounds = rep(Inf, d))