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

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

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.

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.

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
#> $prob
#> [1] 4.180292e-27
#> 
#> $relErr
#> [1] 0.001089263
#> 
#> $upbnd
#> [1] 7.8651e-27
#>