condTruncMVN: Conditional Truncated Multivariate Normal Distribution

Paul M. Hargarten

2020-09-10

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 [1]. 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
#>        [,1]  [,2] [,3]  [,4]   [,5]
#> [1,] 1.0000 0.900 0.81 0.729 0.6561
#> [2,] 0.9000 1.000 0.90 0.810 0.7290
#> [3,] 0.8100 0.900 1.00 0.900 0.8100
#> [4,] 0.7290 0.810 0.90 1.000 0.9000
#> [5,] 0.6561 0.729 0.81 0.900 1.0000

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))
#> $condMean
#> [1]  0.3430923 -0.3211143 -0.8000000
#> 
#> $condVar
#>              [,1]         [,2] [,3]
#> [1,] 1.394510e-01 6.934025e-02 0.00
#> [2,] 6.934025e-02 1.394510e-01 0.00
#> [3,] 1.110223e-16 1.110223e-16 0.19
#> 
#> $condLower
#> [1] -10 -10 -10
#> 
#> $condUpper
#> [1] 10 10 10
#> 
#> $condInit
#> [1] 0 0 0

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)
#> [1] -3.07231

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))
#> [1] 0.01306111
#> attr(,"error")
#> [1] 7.108294e-07
#> attr(,"msg")
#> [1] "Normal Completion"

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))
#>            [,1]       [,2]       [,3]
#> [1,] 0.02238382 -0.1426882 -1.7914223
#> [2,] 0.32684386 -0.5239659 -0.1189072

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))
#> univariate CDF: using truncnorm::ptruncnorm
#> [1] 7.080093e-08

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))
#> [1] 1.160074 1.040832

Computational Details

This vignette is successfully processed using the following.

#>  -- Session info ---------------------------------------------------
#>  setting  value                       
#>  version  R version 4.0.2 (2020-06-22)
#>  os       macOS High Sierra 10.13.6   
#>  system   x86_64, darwin17.0          
#>  ui       X11                         
#>  language (EN)                        
#>  collate  C                           
#>  ctype    en_US.UTF-8                 
#>  tz       America/New_York            
#>  date     2020-09-10
#> --  Packages -------------------------------------------------------
#>  package      * version date       lib source                                 
#>  condMVNorm     2020.1  2020-03-18 [3] CRAN (R 4.0.2)                         
#>  matrixNormal   0.0.4   2020-08-26 [3] CRAN (R 4.0.2)                         
#>  tmvmixnorm     1.1.0   2020-05-19 [3] CRAN (R 4.0.2)                         
#>  tmvtnorm       1.4-10  2015-08-28 [3] CRAN (R 4.0.2)                         
#>  truncnorm      1.0-8   2020-07-27 [3] Github (olafmersmann/[email protected])
#> 
#> [1] /private/var/folders/55/n1q58r751yd7x4vqzdc8zk_w0000gn/T/RtmpG3NqEi/Rinst11cfe5e82216d
#> [2] /private/var/folders/55/n1q58r751yd7x4vqzdc8zk_w0000gn/T/RtmpWVLW4l/temp_libpath1170138ff3692
#> [3] /Library/Frameworks/R.framework/Versions/4.0/Resources/library

Final Notes

Lifecycle::experimental CRAN status

References

1. Horrace, W.C. Some results on the multivariate truncated normal distribution. Journal of Multivariate Analysis 2005, 94, 209–221.