markovianZigzag: Markovian Zigzag Sampler

Description

Sample from a truncated multivariate normal distribution using the Markovian Zigzag process, a continuous-time, non-reversible Markov chain Monte Carlo method based on piecewise deterministic Markov processes (PDMPs).

Usage

markovianZigzag(
  nSample,
  burnin = 0,
  mean,
  prec,
  lowerBounds,
  upperBounds,
  init = NULL,
  stepSize = NULL,
  seed = 1,
  diagnosticMode = FALSE,
  nStatusUpdate = 0L
)

Value

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

Arguments

nSample: Number of samples after burn-in.
burnin: Number of burn-in samples (default = 0).
mean: A d-dimensional mean vector.
prec: A d-by-d precision matrix of the Gaussian distribution.
lowerBounds: A d-dimensional vector specifying the lower bounds. -Inf is accepted.
upperBounds: A d-dimensional vector specifying the upper bounds. Inf is accepted.
init: A d-dimensional vector of the initial value. init must satisfy all constraints. If init = NULL, a random initial value will be used.
stepSize: Step size for the Markovian Zigzag sampler. Default value is the empirically optimal choice: \(\sqrt{2}\lambda^{-1/2}\), where \(\lambda\) is the minimal eigenvalue of the precision matrix.
seed: Random seed (default = 1).
diagnosticMode: Logical. TRUE for also returning diagnostic information such as the step size used.
nStatusUpdate: Number of status updates to print during sampling. If 0 (default), no updates are printed.

References

Bierkens, J., Roberts, G. O., and Zitt, P.-A. (2019). Ergodicity of the zigzag process. The Annals of Applied Probability, 29(4): 2266-2301.

Examples

Run this code

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