recintab: Moments of arbitrary order of a (possibly) truncated multivariate normal variable

Description

Produces an array with the moments up to specified orders of a (possibly) truncated multivariate normal distribution. Each component variable can be truncated on one side (to the left or to the right) or on two sides or not truncated.

Usage

recintab(kappa, a, b, mu, S, ...)

Arguments

kappa

a vector of non-negative integer values representing the required order of moments for each component variable.

a vector representing the lower truncation values of the component variables; -Inf values are allowed.

a vector representing the upper truncation values of the component variables; Inf values are allowed.

a vector representing the mean value of the pre-truncation normal random variable.

a symmetric positive-definite matrix representing the variance matrix of the pre-truncation normal random variable.

…

parameters passed to sadmvn; see the ‘Details’.

Value

In the multivariate case, for an input vector kappa=c(k1,…, kd), the functions returns an array of dimension c((k1+1),...,(kd+1)) whose entries represent integrals described in section ‘Details’. In other words, the array element M[i+1, j+1, k+1,...] contains the unnormalized cross moment of order (i, j, k,...); this must be divided by M[1] to obtain the regular cross moment.

In the univariate case, a vector is returned with similar meaning.

Warning

Although the underlying algorithm is exact in principle, the actual computation hinges crucially on the initial integration of the multivariate normal density over the truncation hyper-cube. This integration may result in numerical inaccuracies, whose amount depends on the supplied arguments. Moreover, the recursion employed by the algorithm propagates the initial error to other terms.

When problematic cases have been processed by the original Matlab function, the same issues have occurred, up to minor variations.

Instances of such errors may be detected when the array M/M[1] is passed to mom2cum, but there is no guarantee that all such problems are detected.

Details

This function is the R translation of the Matlab function with the same name belonging to the package ftnorm, which is associated to the paper of Kan and Robotti (2017). The Matlab package ftnorm has been downloaded from http://www-2.rotman.utoronto.ca/~kan/research.htm, on 2020-04-23.

The function returns an array, M say, whose entries represent integrals of type \(\int_a^b x^\kappa f(x) dx\), where \(f(x)\) denotes the \(d\)-dimensional normal density function. Typically, interest is in the transformed array M/M[1] whose entries represent moments of the truncated distribution.

The algorithm is based on a recursion starting from the integral of the normal distribution over the specified hyper-rectangle. This integral is evaluated by sadmvn, whose tuning parameters maxpts, abseps, releps can be regulated via the ... argument.

References

Kan, Raymond and Robotti, Cesare (2017). On moments of folded and truncated multivariate normal distributions. Journal of Computational and Graphical Statistics, 26, 930-934, DOI: 10.1080/10618600.2017.1322092

Leppard, P. and Tallis, G. M. (1989). Algorithm AS249: Evaluation of the mean and covariance of the truncated multinormal distribution Applied Statistics 38, 543-553)

Examples

Run this code

# NOT RUN {
mu <- c(1, -0.5, 0)
Sigma <- toeplitz(1/(1:3))
low <- c(-Inf, -3, -4)
hi <- c(1.5, Inf, 2)
M <- recintab(c(2,3,1), low, hi, mu, Sigma)
M/M[1]
# cross-moments up to order 2 for X1, up to the 3 for X2, up to 1 for X3,
# if the components of the trivariate variable are denoted (X1,X2,X3)
#--
# Example 2 of Leppard & Tallis (1989, Appl.Stat. vol.38, p.547)
truncp <- c(0, 1, 2)
U <- c(0, 0, 0)
V <- 0.5*(diag(3) + matrix(1, 3, 3))
M <- recintab(c(2,2,2), truncp, rep(Inf,3), U, V)
mom <- M/M[1]
EX <-  c(mom[2,1,1], mom[1,2,1], mom[1,1,2])
print(EX, digits=9)
EX2 <- matrix(c(
          mom[3,1,1], mom[2,2,1], mom[2,1,2],
          mom[2,2,1], mom[1,3,1], mom[1,2,2],
          mom[2,1,2], mom[1,2,2], mom[1,1,3]), 
          3, 3, byrow=TRUE)
varX <- EX2 - outer(EX ,EX)       
print(varX, digits=9)


# }

Run the code above in your browser using DataCamp Workspace