condTruncMVN: Conditional Truncated Multivariate Normal Distribution
In condTruncMVN: Conditional Truncated Multivariate Normal Distribution

knitr::opts_chunk$set(
  collapse = TRUE,  # If TRUE, all output would be in the code chunk.
  comment = "#>",
  fig.path = "man/figures/README-",
  out.width = "100%",
  results = "markup",
  prompt = TRUE,
  strip.white = TRUE,
  tidy = TRUE,
  tidy.opts = list(width.cutoff = 90), # options for tidy to remove blank lines [blank = FALSE] and set the approximate line width to be 80.
  fig.show = "asis",
  fig.height = 4.5,  # inches
  fig.width = 4.5   # inches
)
library("formatR")

The goal of condTruncMVN is to find densities, probabilities, and samples from a conditional truncated multivariate normal distribution. Suppose that Z = (X,Y) is from a fully-joint multivariate normal distribution of dimension n with mean $\boldsymbol\mu$ and covariance matrix sigma ( $\Sigma$ ) truncated between lower and upper. Then, Z has the density $$ f_Z(\textbf{z}, \boldsymbol\mu, \Sigma, \textbf{lower}, \textbf{upper})= \frac{exp(-\frac{1}{2}(\textbf{z}-\boldsymbol\mu)^T \Sigma^{-1} (\textbf{z}-\boldsymbol\mu))} {\int_{\textbf{lower}}^{\textbf{upper}} exp(-\frac{1}{2}(\textbf{z}-\boldsymbol\mu)^T \Sigma^{-1} (\textbf{z}-\boldsymbol\mu)) d\textbf{z}} $$ for all z in [lower, upper] in $\mathbb{R^{n}}$.

This package computes the conditional truncated multivariate normal distribution of Y|X. The conditional distribution follows a truncated multivariate normal [@horraceResultsMultivariateTruncated2005]. Specifically, the functions are arranged such that

$$ Y = Z[ , dependent.ind] $$ $$ X = Z[ , given.ind] $$ $$ Y|X = X.given \sim MVN(mean, sigma, lower, upper) $$ The [d,p,r]cmvtnorm() functions create a list of parameters used in truncated conditional normal and then passes the parameters to the source function below.

| Function Name | Description | Source Function | Univariate Case| Additional Parameters |---------- | ------------| ---------- | ------------ | --------------------- |
| condtMVN | List of parameters used in truncated conditional normal.| condMVNorm:: condMVN() | | | | dcmvtnorm | Calculates the density f(Y=y\| X = X.given) up to a constant. The integral of truncated distribution is not computed. | tmvtnorm:: dtmvnorm() | truncnorm:: dtruncnorm() | y, log | | pcmvtnorm | Calculates the probability that Y\|X is between lowerY and upperY given the parameters. | tmvtnorm:: ptmvnorm() | truncnorm:: ptruncnorm() | lowerY, upperY, maxpts, abseps, releps | | rcmvtnorm | Generate random sample. | tmvmixnorm:: rtmvn() | truncnorm:: rtruncnorm() | n, init, burn, thin |

Installation

You can install the released version of condTruncMVN from CRAN with install.packages("condTruncMVN"). You can load the package by:

library("condTruncMVN")

And the development version from GitHub with:

Example

Suppose $X2,X3,X5|X1 = 1,X4 = -1 \sim N_3(1, Sigma, -10, 10)$. The following code finds the parameters of the distribution, calculates the density, probability, and finds random variates from this distribution.

library(condTruncMVN)
d <- 5
rho <- 0.9
Sigma <- matrix(0, nrow = d, ncol = d)
Sigma <- rho^abs(row(Sigma) - col(Sigma))
Sigma

First, we find the conditional Truncated Normal Parameters.

condtMVN(mean  = rep(1, d),
  sigma = Sigma,
  lower = rep(-10, d),
  upper = rep(10, d),
  dependent.ind = c(2, 3, 5),
  given.ind = c(1, 4), X.given = c(1, -1)
)

Find the log-density when X2,X3,X5 all equal $0$:

dcmvtruncnorm(
  rep(0, 3),
  mean  = rep(1, 5),
  sigma = Sigma,
  lower = rep(-10, 5),
  upper = rep(10, d),
  dependent.ind = c(2, 3, 5),
  given.ind = c(1, 4), X.given = c(1, -1),
  log = TRUE
)

Find $P( -0.5 < X2,X3,X5 < 0 | X1 = 1,X4 = -1)$:

pcmvtruncnorm(
  rep(-0.5, 3), rep(0, 3),
  mean = rep(1, d),
  sigma = Sigma,
  lower = rep(-10, d),
  upper = rep(10, d),
  dependent.ind = c(2, 3, 5),
  given.ind = c(1, 4), X.given = c(1, -1)
)

Generate two random numbers from the distribution.

set.seed(2342)
rcmvtruncnorm(2,
  mean = rep(1, d),
  sigma = Sigma,
  lower = rep(-10, d),
  upper = rep(10, d),
  dependent.ind = c(2, 3, 5),
  given.ind = c(1, 4), X.given = c(1, -1)
)

Another Example: To find the probability that $X1|X2, X3, X4, X5 \sim N(1, Sigma, -10, 10)$ is between -0.5 and 0:

pcmvtruncnorm(-0.5, 0,
  mean  = rep(1, d),
  sigma = Sigma,
  lower = rep(-10, d),
  upper = rep(10, d),
  dependent.ind = 1,
  given.ind = 2:5, X.given = c(1, -1, 1, -1)
)

If I want to generate 2 random variates from $X1|X2, X3, X4, X5 \sim N(1, Sigma, -10, 10)$:

set.seed(2342)
rcmvtruncnorm(2,
  mean  = rep(1, d),
  sigma = Sigma,
  lower = rep(-10, d),
  upper = rep(10, d),
  dependent.ind = 1,
  given.ind = 2:5, X.given = c(1, -1, 1, -1)
)

Computational Details

This vignette is successfully processed using the following.

# sessioninfo::session_info() # makes a mess! Instead

cat(" -- Session info ---------------------------------------------------")
sessioninfo::platform_info()
cat("--  Packages -------------------------------------------------------")
tmp.df <- sessioninfo::package_info(
  c("condMVNorm", "matrixNormal", "tmvmixnorm", "tmvtnorm", "truncnorm"),
  dependencies = FALSE
)
print(tmp.df)