computes an estimator and a deterministic upper bound of the probability Pr(l<X<u), where X is a zeromean multivariate normal vector with covariance matrix Σ, that is, X is drawn from N(0,Σ) infinite values for vectors u and l are accepted; Monte Carlo method uses sample size n;
1  mvNcdf(l, u, Sig, n)

l 
lower truncation limit 
u 
upper truncation limit 
Sig 
covariance matrix of N(0,Σ) 
n 
Monte Carlo simulation effort — the larger the n, the smaller the relative error of the estimator. 
Suppose you wish to estimate p=Pr(l<AX<u), where A is a full rank matrix and X is drawn from N(μ,Σ), then you simply compute p=Pr(lAμ<AY<uAμ), where Y is drawn from N(0, AΣ A^\top).
est
with structure

estimated value of probability Pr(l<X<u) 

estimated relative error of estimator 

theoretical upper bound on true Pr(l<X<u) 
For small dimensions, say d<50, better accuracy may be obtained by using
the (usually slower) quasiMonte Carlo version mvNqmc
of this algorithm.
Z. I. Botev, email: botev@unsw.edu.au and web page: http://web.maths.unsw.edu.au/~zdravkobotev/
Z. I. Botev (2015), The Normal Law Under Linear Restrictions: Simulation and Estimation via Minimax Tilting, submitted to JRSS(B)
See also: mvNqmc
and mvrandn
1 2 3 4 5 
Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.
All documentation is copyright its authors; we didn't write any of that.