mvNqmc: Truncated multivariate normal cumulative distribution...

In lbelzile/TruncatedNormal: Truncated Multivariate Normal and Student Distributions

Truncated multivariate normal cumulative distribution (quasi-Monte Carlo)

Description

Computes an estimate and a deterministic upper bound of the probability Pr(l<X<u), where X is a zero-mean multivariate normal vector with covariance matrix \Sigma, that is, X is drawn from N(0,\Sigma). Infinite values for vectors u and l are accepted. The Monte Carlo method uses sample size n: the larger n, the smaller the relative error of the estimator.

Usage

mvNqmc(l, u, Sig, n = 1e+05)

Arguments

l	lower truncation limit
u	upper truncation limit
Sig	covariance matrix of N(0,\Sigma)
n	number of Monte Carlo simulations

Details

Suppose you wish to estimate Pr(l<AX<u), where A is a full rank matrix and X is drawn from N(\mu,\Sigma), then you simply compute Pr(l-A\mu<AY<u-A\mu), where Y is drawn from N(0, A\Sigma A^\top).

Value

a list with components

• prob: estimated value of probability Pr(l<X<u)

• relErr: estimated relative error of estimator

• upbnd: theoretical upper bound on true Pr(l<X<u)

Note

This version uses a Quasi Monte Carlo (QMC) pointset of size ceiling(n/12) and estimates the relative error using 12 independent randomized QMC estimators. QMC is slower than ordinary Monte Carlo, but is also likely to be more accurate when d<50. For high dimensions, say d>50, you may obtain the same accuracy using the (typically faster) mvNcdf.

Author(s)

Zdravko I. Botev

References

Z. I. Botev (2017), The Normal Law Under Linear Restrictions: Simulation and Estimation via Minimax Tilting, Journal of the Royal Statistical Society, Series B, 79 (1), pp. 1–24.

See Also

mvNcdf, mvrandn

Examples

d <- 15 
l <- 1:d
u <- rep(Inf, d)
Sig <- matrix(rnorm(d^2), d, d)*2 
Sig <- Sig %*% t(Sig)
mvNqmc(l, u, Sig, 1e4) # compute the probability

lbelzile/TruncatedNormal documentation built on March 4, 2024, 5:50 p.m.