Problem with mvncdf - returning NaN

Hello, I am having some problems when using 'mvncdf'. Consider:

mu = [0.4704 14.5809];

cov = [0.0030 0.0922; 0.0922 2.8597];

xl = [8 0];

xu = [12 4];

Then, mvncdf(xl, xu, mu, cov) returns NaN but I cannot understand why. What is more, the documentation doesn't mention anything relative. And, what is even more strange, mvncdf(xu, mu, cov) - mvncdf(xl, mu, cov) (which is equivalent) returns a double, instead of Nan. Is that due to a bug in mvncdf, or do I miss something important?

Answer 1

[Paul]

The limits of area of integration defined in xl and xu are way, way, out in the tails of the density function, which may be why mvncdf is having an issue.

As I understand mvncdf ...let p(x1,x2) be the joint density. Then mvncdf computes the probability, P, given by:

syms p(x1,x2)
xl = sym('xl',[1 2]);
xu = sym('xu',[1 2]);
P = int(int(p,x1,xl(1),xu(1)),x2,xl(2),xu(2))

P = 

First, let's take a look at the pdf as defined by the parameters in the question

mu = [0.4704   14.5809];
cov = [0.0030    0.0922;    0.0922    2.8597];
func(mu,cov);

Now, redraw but with a rectangle showing the area of integration:

xlower = [8 0];  
xupper = [12 4];
func(mu,cov,xlower,xupper)

As you can see, the area of integration is way out there. I don't know why mvncdf returns NaN in this case instead of zero, but it probably has something to do with trying to integrate a function with incredibly small values.

Now try with a more reasonable area of integration

xlower = [0.4 13];
xupper = [0.5 15];
func(mu,cov,xlower,xupper);

The probability is:

mvncdf(xlower,xupper,mu,cov)

ans = 0.4229

which seems pretty reasoable.

Also, I don't think that mvncdf(xu, mu, cov) - mvncdf(xl, mu, cov) is equivalent

mvncdf(xupper, mu, cov) - mvncdf(xlower, mu, cov)

ans = 0.4985

because it includes area outside of the red rectangle, particularly the light blue region where x1 > 0.4 and x2 < 13.

function func(mu,cov,xlower,xupper)
x1 = 0.2:.001:.8;
x2 = 10:.001:18;
[X1,X2] = meshgrid(x1,x2);
X = [X1(:) X2(:)];
p = mvnpdf(X,mu,cov);
p = reshape(p,length(x2),length(x1));
figure
pcolor(x1,x2,p),shading interp,colorbar
xlabel('x1'),ylabel('x2')
if nargin == 4
rectangle('Position',[xlower(1) xlower(2) xupper(1)-xlower(1) xupper(2)-xlower(2)],'EdgeColor','r');
axis auto
end
end

[Georgios Apostolakis]

Oh I see, thank you very much!

Thus, is it safe to replace all NaN returned by 'mvncdf (xupper, xlower, mu, cov)' with 0?

[Paul]

I don't know if that's safe or not, because I don't know exactly why mvncdf is returning NaN. Depending on what you're doing, it may be safe for you to make a judgement about when to assign the result to zero based on the area of integration and how far away it is (relative to cov) from the mean.

It might also be worthwhile to reach out to Tech Support about this, because mvncdf shouldn't be returning NaN.

[Georgios Apostolakis]

Tech support of Matlab is very unreliable. I cannot reach out to them as I use R2014a for which there is no support anymore.

Despite the existence of the same problem in R2021a, they refuse to check it out (and fix it in future versions) unless I provide them a valid license number of a newer Matlab version (which obviously I have not, the only valid license I own is the one of R2014a, provided by my university).

If anyone reading this post owns a valid license of a relatively new Matlab version, let me know so I can report this bug to tech support. Or, he can report it himself if he prefers. Thank you!

[Paul]

Submitted. I'll repport back here when I get a response.

[Georgios Apostolakis]

Thanks!!

[Paul]

I submitted against 2020b. The response was: "Our developer team confirmed the issue, and they will consider a resolution for a future release of MATLAB."

[Georgios Apostolakis]

Great. Thank you!!

[adam]

FWIW, this bug is still present in my current version of MATLAB 2021a:

>> x

x =

0.7551 0.4435 -1.5760

>> mu

mu =

-0.8469 -1.9076 2.3800

>> covar

covar =

0.0052 0.0045 0.0027

0.0045 0.0063 0.0028

0.0027 0.0028 0.0015

>> mvncdf(x,mu,covar)

ans =

NaN

[Paul]

Fixed in either 2021b or 2022a, it appears, unless the extra digits not shown make a difference. However, I didn't get a notice of such, which I have for other issues that I've reported that have been fixed.

x = [0.7551    0.4435   -1.5760];
mu = [-0.8469   -1.9076    2.3800];
covar = [
    0.0052    0.0045    0.0027
    0.0045    0.0063    0.0028
    0.0027    0.0028    0.0015];

mvncdf(x,mu,covar)

ans = 0