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.
recintab(kappa, a, b, mu, S, ...)
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
a vector representing the upper truncation values of the component
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
The maximal dimension of the multivariate normal variable is 20.
If this threshold is exceeded
NAs are returned.
This function is the R translation of the Matlab function with the same name
belonging to the package
which is associated to the paper of Kan and Robotti (2017).
The Matlab package
ftnorm has been downloaded from
The function returns an array,
M say, whose entries represent integrals
of type \int_a^b x^κ f(x) dx, where f(x) denotes the
d-dimensional normal density function.
Typically, interest is in the scaled array
M/M whose entries
represent the 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
sadmvn, whose tuning parameters
maxpts, abseps, releps
can be regulated via the
In the multivariate case, for an input vector
the functions returns an array of dimension
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 to obtain the regular cross moment.
In the univariate case, a vector is returned with similar meaning.
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
mom2cum, but there is no guarantee that all such
problems are detected.
This function is not intended for direct call by a user, at least in
commonly encountered situations.
mom.mtruncnorm represents a more user-friendly tool.
Original Matlab code by Raymond Kan and Cesare Robotti, porting to R by Adelchi Azzalini.
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)
mom.mtruncnorm for a more user-friendly function,
mom2cum for transformation to cumulants,
sadmvn for regulating accuracy if
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 # 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 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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.