mcmcsae (version 0.6.0)

create_TMVN_sampler: Set up a sampler object for sampling from a possibly truncated and degenerate multivariate normal distribution

Description

This function sets up an object for multivariate normal sampling based on a specified precision matrix or cholesky decomposition thereof. Linear equality and inequality restrictions are supported. For sampling under inequality restrictions three algorithms are available. The default in that case is an exact Hamiltonian Monte Carlo algorithm (Pakman and Paninski, 2014). Alternatively, a Gibbs sampling algorithm can be used (Rodriguez-Yam et al., 2004). The third option is a data augmentation method to sample from a smooth approximation to the truncated multivariate normal distribution (Souris et al., 2018).

Usage

create_TMVN_sampler(
  Q,
  perm = NULL,
  mu = NULL,
  Xy = NULL,
  update.Q = FALSE,
  update.mu = update.Q,
  name = "x",
  coef.names = NULL,
  R = NULL,
  r = NULL,
  S = NULL,
  s = NULL,
  lower = NULL,
  upper = NULL,
  method = NULL,
  reduce = (method == "Gibbs" && !is.null(R)),
  T.HMC = pi/2,
  print.info = FALSE,
  sharpness = 100,
  useV = FALSE
)

Arguments

Q

precision matrix (unconstrained) (n x n); used in case cholQ is not supplied.

perm

whether permutation/pivoting should be used to build a Cholesky object.

mu

mean of the (unconstrained) multivariate normal distribution.

Xy

alternative to specifying mu; in this case mu is computed as Q^{-1} Xy.

update.Q

whether Q is updated for each draw.

update.mu

whether mu is updated for each draw. By default equal to update.Q.

name

name of the TMVN vector parameter.

coef.names

optional labels for the components of the vector parameter.

R

equality restriction matrix.

r

rhs vector for equality constraints R^T x = r.

S

inequality restriction matrix.

s

rhs vector for inequality constraints S^T x >= s.

lower

alternative to s for two-sided inequality restrictions lower <= S^T x <= upper.

upper

alternative to s for two-sided inequality restrictions lower <= S^T x <= upper.

method

sampling method. The options are "direct" for direct sampling from the unconstrained or equality constrained multivariate normal (MVN). For inequality constrained MVN sampling three methods are supported: "HMC" for (exact) Hamiltonian Monte Carlo, "Gibbs" for a component-wise Gibbs sampling approach, and "softTMVN" for a data augmentation method that samples from a smooth approximation to the truncated MVN.

reduce

whether to a priori restrict the simulation to the subspace defined by the equality constraints.

T.HMC

the duration of a Hamiltonian Monte Carlo simulated particle trajectory. If a vector of length 2 is supplied it is interpreted as an interval from which the duration is drawn uniformly, independently in each HMC iteration.

print.info

whether information about violations of inequalities and bounces off inequality walls is printed to the screen. This sometimes provides useful diagnostic information.

sharpness

for method 'softTMVN', the sharpness of the soft inequalities; the larger the better the approximation of exact inequalities. It must be either a scalar value or a vector of length equal to the number of inequality restrictions.

useV

for method 'softTMVN' whether to base computations on variance instead of precision matrices.

Value

An environment for sampling from a possibly degenerate and truncated multivariate normal distribution.

Details

The componentwise Gibbs sampler uses univariate truncated normal samplers as described in Botev and L'Ecuyer (2016). These samplers are implemented in R package TruncatedNormal, but here translated to C++ for an additional speed-up.

References

Z.I. Botev and P. L'Ecuyer (2016). Simulation from the Normal Distribution Truncated to an Interval in the Tail. in VALUETOOLS.

Y. Cong, B. Chen and M. Zhou (2017). Fast simulation of hyperplane-truncated multivariate normal distributions. Bayesian Analysis 12(4), 1017-1037.

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

G. Rodriguez-Yam, R.A. Davis and L.L. Scharf (2004). Efficient Gibbs sampling of truncated multivariate normal with application to constrained linear regression. Unpublished manuscript.

H. Rue and L. Held (2005). Gaussian Markov Random Fields. Chapman & Hall/CRC.

A. Souris, A. Bhattacharya and P. Debdeep (2018). The Soft Multivariate Truncated Normal Distribution. arXiv:1807.09155.

Examples

# NOT RUN {
S <- cbind(diag(2), c(-1, 1), c(1.1, -1))  # inequality matrix
# S'x >= 0 represents the wedge x1 <= x2 <= 1.1 x1
# example taken from Pakman and Paninski (2014)
# 1. exact Hamiltonian Monte Carlo (Pakman and Paninski, 2014)
sampler <- create_TMVN_sampler(Q=diag(2), mu=c(4, 4), S=S, method="HMC")
sim <- MCMCsim(sampler)
summary(sim)
plot(as.matrix(sim$x), pch=".")
# 2. Gibbs sampling approach (Rodriguez-Yam et al., 2004)
sampler <- create_TMVN_sampler(Q=diag(2), mu=c(4, 4), S=S, method="Gibbs")
sim <- MCMCsim(sampler)
summary(sim)
plot(as.matrix(sim$x), pch=".")
# 3. soft TMVN approximation (Souris et al., 2018)
sampler <- create_TMVN_sampler(Q=diag(2), mu=c(4, 4), S=S, method="softTMVN")
sim <- MCMCsim(sampler)
summary(sim)
plot(as.matrix(sim$x), pch=".")
# }
# NOT RUN {
# }