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 |
You can install the released version of condTruncMVN from CRAN with install.packages("condTruncMVN")
. You can load the package by:
And the development version from GitHub with:
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**)\):
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/truncnorm@eea186e)
#>
#> [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
1. Horrace, W.C. Some results on the multivariate truncated normal distribution. Journal of Multivariate Analysis 2005, 94, 209–221.