Package 'trunmnt'

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, [email protected]

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

Obj	An mtrunmnt object created by the mtrunmnt.

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

mu	**(Required)** Mean vector of the parent multivariate normal distribution.
lower	vector of lower limits. If the lower limits are the same, a scalar value can be given. Defaults to -Inf.
upper	Vector of upper limits. If the upper limits are the same, a scalar value can be given. Defaults to Inf.
Sigma	The variance-covariance matrix of the parent multivariate normal distribution. It must be given a symmetric positive definite matrix, if Sigmae, D and Z are not specified.
Sigmae	Vector of variances of error terms. If the variances are the same, a scalar value can be given. Defaults to 1.
Z	Design matrix for the random components. Defaults to $n \times 1$ matrix of 1's where $n$ is the dimension of mu. It must be specified carefully with the argument D. ncol(Z) determines the dimension of D.
D	Variance-covariance matrix of $u$. See Details. If the random components are independent, you can specify either a vector of variances or a scalar value. A scalar value means that the random components have the same variance. Defaults to 1.
nGH	Number of quadrature points. Defaults to 35.

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

Obj	An mtrunmnt object created by the mtrunmnt.

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

Obj	mtrunmnt object created by the mtrunmnt.
kappa	Vector of orders of length equal to mu.

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

k	Order of moment. It must be a non-negative integer.
mu	Mean of parent normal distribution. Defaults to 0.
lower	Lower limit. Defaults to -Inf.
upper	Upper limit. Defaults to Inf.
sd	Standard deviation of parent normal distribution. Defaults to 1.

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)

---