mvNqmc function

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 $\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.

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

Returns

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)$

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)$.

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.

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

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

Author(s)

Zdravko I. Botev