trunmnt (version 1.0.0)

trunmnt-package: 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.

Author

Seung-Chun Lee,

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.

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)