File Exchange

image thumbnail

Multivariate normal cumulative distribution

version 1.1.0.0 (9.9 KB) by Zdravko Botev
state-of-the-art algorithm for computing the multivariate normal cdf in high dimensions

4 Downloads

Updated 28 Oct 2015

View Version History

View License

Computes the probability Pr(l<X<u), where 'X' is a zero-mean multivariate normal vector with covariance 'Sig'.
In high dimensions, this algorithm is vastly superior to the one in Matlab's statistics toolbox, see example.
Reference: Z. I. Botev (2015), "The Normal Law Under Linear Restrictions: Simulation and Estimation via Minimax Tilting", submitted to JRSS(B)

Cite As

Zdravko Botev (2020). Multivariate normal cumulative distribution (https://www.mathworks.com/matlabcentral/fileexchange/53583-multivariate-normal-cumulative-distribution), MATLAB Central File Exchange. Retrieved .

Comments and Ratings (9)

Lu Xiaohui

very useful, Solve my problems when d>25. Thanks!

Qi Wu

man, I should rate you 5-star even before I use the code. This 25-limitation of the inbuild function in matlab is useless

LuisCardona

Zdravko Botev

Another example: in the code below, if you change to
Sig=inv(0.5*eye(d)+.5*ones(d,d));
then my code is more than 10,000 times more efficient!
I hope this performance will earn a 5-star rating ;-)

Zdravko Botev

Hi LuisCardona,
yet again, you are not using the code appropriately, and are making incorrect conclusions.
Please read the documentation carefully if you want to obtain better efficiency.

I do not know how much you know about the difference between quasi Monte Carlo and Monte Carlo.
For the dimensions you have chosen, you have to use Quasi Monte Carlo (mvNqmc), not Monte Carlo
(mvNcdf). Quasi Monte Carlo is more efficient in low dimensions, but becomes uncompetitive in high-dimensions (Matlab’s function does not go beyond d=25, whereas my code can handle d>25)

If you use the code correctly, then my code is about 20 times more efficient, as demonstrated here:

clear all,clc,d=9;
l=ones(d,1)/2;u=ones(d,1);Sig=(0.5*eye(d)+.5*ones(d,d));
tic
for i=1:15
options=optimset('TolFun',0,'MaxFunEvals',10^6);
prob(i)=mvncdf(l,u,zeros(d,1),Sig,options);
end
efficiency1=var(prob)*toc
tic
for i=1:15
est=mvNqmc(l,u,Sig,10^6);
e(i)=est.prob;
end
efficiency2=var(e)*toc
efficiency1/efficiency2

LuisCardona

Thank you for clarifying that. I made another test with 9 variables. Like this
tic; for i=1:50; t=mvncdf(l,u,rho9, 5*10^5); est(i)=t.prob; end; toc

The matlab built in took 21.903402 seconds and the standard deviation of the values was 2.4996e-05. This function took 23.656228 seconds and the standard deviation of the values was 1.7669e-04. So, the matlab built had a lower variance.

Zdravko Botev

LuisCardona, the reason you get this result is that you used dimension=d=3, which is too low and does not need Monte Carlo simulation. For d<=3, Matlab uses deterministic numerical quadrature, so naturally there is no randomness and the variance is 0. If you try, say d=5, or higher, then mvNcdf performs much better, because there is no deterministic quadrature that can compute accurately integrals of high dimension.

LuisCardona

The explanation for the results that I obtained in the last comment may be because I = -Inf. May be for finite intervals, it outperforms

LuisCardona

I tested this function calculated 50 times with the same parameters to compare this one and the one from Matlab for 3 random variables.

Like this
for i=1:50; est = mvncdf(l,u,Sig,10^6); devs(i)= est.prob; end; toc

this one took 12.184780 seconds and matlab took 0.039231 seconds. The std of this one was 8.7441e-05, while matlab's was 0.

MATLAB Release Compatibility
Created with R2015b
Compatible with any release
Platform Compatibility
Windows macOS Linux

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!