Truncated multivariate normal cumulative distribution (quasi-Monte Carlo)
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 Σ, that is, X is drawn from N(0,Σ). 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.
mvNqmc(l, u, Sig, n = 1e+05)
l |
lower truncation limit |
u |
upper truncation limit |
Sig |
covariance matrix of N(0,Σ) |
n |
number of Monte Carlo simulations |
Suppose you wish to estimate Pr(l<AX<u), where A is a full rank matrix and X is drawn from N(μ,Σ), then you simply compute Pr(l-Aμ<AY<u-Aμ), where Y is drawn from N(0, AΣ A^\top).
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)
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.
Zdravko I. Botev
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.
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
Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.