Moments of Truncated Multivariate Normal Distribution.

Description

Based on the algorithm given by Lee (2021), it computes the product moment of a truncated multivariate normal distribution using the multivariate Gaussian quadrature.

Details

Use mtrunmnt to generate a S3 objective and then use prodmnt to compute arbitrary order of moments. It can also calculate the mean and variance-covariance matrix of a truncated multivariate normal distribution by meanvar. See also probntrun and utrunmnt for computing the probability of multivariate normal distribution for given limits, and the moment of truncated univariate normal distribution.

The 'trunmnt' package computes the product moments of the Truncated multivariate normal distribution by implementing the algorithm proposed by Lee (2020). This approach relies on multivariate Gaussian quadrature for numerical integration.

**Limitations:** While the package supports arbitrary positive-definite matrices, the computational complexity scales poorly with the dimension of the vector. For dimensions $D > 5$, using the pattern-specific functions or considering approximation methods is highly recommended.

Author(s)

Seung-Chun Lee, seung@hs.ac.kr

References

Burkardt, J. (2014). The truncated normal distribution, Online document, Available from: https://people.sc.fsu.edu/~jburkardt/presentations/truncated_normal.pdf.

Jaeckel, P. (2005). A note on multivariate Gauss-Hermite quadrature. London: ABN-Amro. Avaiable from: http://www.jaeckel.org/ANoteOnMultivariateGaussHermiteQuadrature.pdf.

Lee, S.-C. (2021). Moments Calculation for Truncated Multivariate Normal in Nonlinear Generalized Mixed Models. Communications for Statistical Applications and Methods, 27, 377–383.

See Also

MomTrunc-package

Examples

set.seed(123)
sigma2e <- 1
sigma2a <- 2
n <- 5
mu <- seq(-1,1, length.out = n)
y <- mu + rnorm(1, sd = sqrt(sigma2a)) + rnorm(n, sd = sqrt(sigma2e))
S <- matrix(sigma2a, ncol = n, nrow = n) + diag(sigma2e, n)
a  <- rep(-Inf, n)
b  <- rep(Inf, n)
a[y >= 0] <- 0
b[y <  0] <- 0
obj1 <- mtrunmnt(mu, lower = a, upper = b, Sigmae = sigma2e, D = sigma2a) 
obj2 <- mtrunmnt(mu, lower = a, upper = b, Sigma = S) 
probntrun(obj1)
probntrun(obj2)
prodmnt(obj1, c(2,2,0,0,0))
prodmnt(obj2, c(2,2,0,0,0))
meanvar(obj1)
meanvar(obj2)

Mean and variance for a truncated multivariate normal distribution.

Description

meanvar is a S3 generic function of the class mtrunmnt. Using the prodmnt, it compute the mean and variance-covariance matrix for a truncated multivariate normal distribution.

Usage

meanvar(Obj)

Arguments

Value

A list with the mean and variance-covariance matrix.

See Also

mtrunmnt, meanvarTMD, mtmvnorm and prodmnt.

Examples

###  A simple example ####

set.seed(123)
sigma2e <- 1
sigma2a <- 2
n <- 5
mu <- seq(-1,1, length.out = n)
y <- mu + rnorm(1, sd = sqrt(sigma2a)) + rnorm(n, sd = sqrt(sigma2e))
a  <- rep(-Inf, n)
b  <- rep(Inf, n)
a[y >= 0] <- 0
b[y <  0] <- 0
obj <- mtrunmnt(mu, lower = a, upper = b, Sigmae = sigma2e, D = sigma2a) 
meanvar(obj)

Creating an S3 object for computing the product moment

Description

mtrunmnt creates an S3 object designed to compute the product moment for a truncated multivariate normal distribution, utilizing the algorithm by Lee (2020). Key attributes of this object are the nodes and weights of the multivariate Gaussian quadrature and the probability of the truncation interval

Usage

mtrunmnt(
  mu,
  lower = -Inf,
  upper = Inf,
  Sigma = 1,
  Sigmae = 1,
  Z = matrix(rep(1, length(mu)), ncol = 1),
  D = matrix(1, ncol = 1, nrow = 1),
  nGH = 35
)

Arguments

Details

Assume the parent multivariate normal distribution comes from a mixed-effects linear model:

\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \mathbf{Zu} + \boldsymbol{\epsilon}

where \mathbf{X} and \mathbf{Z} are design matrices corresponding to \boldsymbol{\beta} and \mathbf{u} representing fixed and random effects, respectively, and \boldsymbol{\epsilon} is the vector of errors. It is assumed that the random effects \mathbf{u} follows a multivariate normal distribution with mean \mathbf{0}, and symmetric positive definite variance-covariance matrix \mathbf{D}. As usual, the distribution of \boldsymbol{\epsilon} is assumed to be a multivariate normal with mean \mathbf{0} and variance-covariance matrix \sigma^2_{\epsilon} \mathbf{I}, but for more flexibility, it can be assumed that the error terms are independent, but do not have equal variance. That is, **we assume** \boldsymbol{\epsilon} \sim N(\mathbf{0}, \mathbf{E}) where \mathbf{E} is a diagonal matrix. Then,

\mathbf{Y} \sim N(\boldsymbol{\mu}, \mathbf{\Sigma})

where \mathbf{\Sigma} = \mathbf{Z}\mathbf{D}\mathbf{Z}' + \mathbf{E}. The variance-covariance structure in mtrunmnt can thus be specified either by providing the individual components \mathbf{D}, \mathbf{Z}, and \mathbf{E}, or by directly supplying the resulting overall variance-covariance matrix boldsymbol{\Sigma}.

Value

A mtrunmnt object

References

Lee, S.-C. (2020). Moments calculation for truncated multivariate normal in nonlinear generalized mixed models. Communications for Statistical Applications and Methods, Vol. 27, No. 3, 377–383.

Examples

### Create a mtrunmnt objective ###

set.seed(123)
sigma2e <- 1
sigma2a <- 2
n <- 5
mu <- seq(-1,1, length.out = n)
y <- mu + rnorm(1, sd = sqrt(sigma2a)) + rnorm(n, sd = sqrt(sigma2e))
S <- matrix(sigma2a, ncol = n, nrow = n) + diag(sigma2e, n)
a  <- rep(-Inf, n)
b  <- rep(Inf, n)
a[y >= 0] <- 0
b[y <  0] <- 0
obj1 <- mtrunmnt(mu, lower = a, upper = b, Sigmae = sigma2e, D = sigma2a) 
obj2 <- mtrunmnt(mu, lower = a, upper = b, Sigma = S) 
probntrun(obj1)
probntrun(obj2)
prodmnt(obj1, c(2,2,0,0,0))
prodmnt(obj2, c(2,2,0,0,0))
meanvar(obj1)
meanvar(obj2)

Compute the probability mass after truncation.

Description

Compute the probability for the truncation interval (\mathbf{a}, \mathbf{b}).

Usage

probntrun(Obj)

Arguments

Details

probntrun is a S3 generic function of the class mtrunmnt. Using the multivariate Gaussian quadrature, it computes the probability of truncation interval (\mathbf{a}, \mathbf{b}).

Value

a numeric value.

See Also

mtrunmnt, ptmvnorm.

Examples

###  A simple example ####
set.seed(123)
sigma2e <- 1
sigma2a <- 2
n <- 5
mu <- seq(-1,1, length.out = n)
y <- mu + rnorm(1, sd = sqrt(sigma2a)) + rnorm(n, sd = sqrt(sigma2e))
a  <- rep(-Inf, n)
b  <- rep(Inf, n)
a[y >= 0] <- 0
b[y <  0] <- 0
obj <- mtrunmnt(mu, lower = a, upper = b, Sigmae = sigma2e, D = sigma2a) 
probntrun(obj)

Compute the product moment.

Description

It computes the kappa-th order product moment for a truncated multivariate normal distribution.

Usage

prodmnt(Obj, kappa)

Arguments

Details

prodmnt is a S3 generic function of the class mtrunmnt. Using the multivariate Gaussian quadrature, it computes the product moment E(\prod_{i=1}^n Y_i^{k_i}) | a_i < Y_i < b_i, i = 1, \ldots, n), where Y \sim N(\mu, ZDZ' + E).

Value

a numeric value.

See Also

mtrunmnt, momentsTMD.

Examples

###  A simple example ####

set.seed(123)
sigma2e <- 1
sigma2a <- 2
n <- 5
mu <- seq(-1,1, length.out = n)
y <- mu + rnorm(1, sd = sqrt(sigma2a)) + rnorm(n, sd = sqrt(sigma2e))
a  <- rep(-Inf, n)
b  <- rep(Inf, n)
a[y >= 0] <- 0
b[y <  0] <- 0
obj <- mtrunmnt(mu, lower = a, upper = b, Sigmae = sigma2e, D = sigma2a) 
prodmnt(obj, c(2,2,0,0,0))

Compute the k-th order moment for an univariate truncated normal distribution.

Description

The utrunmnt function uses the moment generating function to compute any order of moment for the truncated normal distribution.

Usage

utrunmnt(k, mu = 0, lower = -Inf, upper = Inf, sd = 1)

Arguments

Value

a numeric value.

References

Burkardt, J. (2014). The truncated normal distribution, Online document, Available from: https://people.sc.fsu.edu/~jburkardt/presentations/truncated_normal.pdf.

Examples

utrunmnt(4, mu = 5, upper = 10)
utrunmnt(1, mu = 5, lower = -3, upper = 4, sd = 2)