Error using mvncdf: "SIGMA must be a square, symmetric, positive definite matrix."

When running mvncdf I get an error: "SIGMA must be a square, symmetric, positive definite matrix." Here is my code:

x = [.125,.125,.125,.125,.125,.125,.125,.195,.195,.11,.135,.135,.145];
z = [466.087,480,483.4783,486.9565,466.087,424.3478,...
463.3043,468.2171,390.6977,335.7143,417.094,478.6325,450.4202];

sigx = std(x);
sigz = std(z);
cov = E((x-E(x)).*(z-E(z)));
cor = cov/(sigx*sigz);

X = [1,1];

mu = [mean(x) mean(z)];
sigma = [sigx^2, cor*sigx*sigz; cor*sigx*sigz, sigz^2];
Y = mvncdf([0,0],X,mu,sigma);

Where E is a short function to find the expected value of a random variable:

function expect = E(q)
expect = 0;
for i=1:length(q)
    expect = expect + q(i)*normpdf(q(i),mean(q),std(q));
end
end

What is causing this error? As far as I can tell everything should be fine with my sigma matrix.

Answer 1

[Walter Roberson]

The underlying math libraries have changed over time, with the exact time of the change depending on the release. For OS-X the change was between R2015b and R2016a. This had an effect on the output of qr() which in turn had an effect on the output of chol(), which is what mvncdf used to test whether the matrix is positive definite. Matrices that were near the boundary of being positive definite might now be calculated as being non positive-definite. You can read a bit more about that at https://www.mathworks.com/matlabcentral/answers/319225-hac-results-vary-bewteen-matlab-r2015b-and-matlab-r2016b

Due to the way that floating point rounding works, you cannot really say that the determination is wrong, even if you can show that if you had used rational numbers to indefinite precision that the matrix was positive definite -- because floating point calculations never do use that.

That said:

I just ran the calculation in symbolic mode; outputting sigma is rather slow, because the expression is pretty long. The expression involves a bunch of exp() of fractions, and all of those are of course irrational numbers... calculating the decomposition is probably not too practical.

To about 50 decimal places your sigma comes out as

[ 0.00071089743589743589743589743589743589743589743589744, -36.493857123342967514465689999208360570489097569933]
[    -36.493857123342967514465689999208360570489097569933,   1952.922555278590113392548958542728033974487558748]

If you explore the form

[x -36; -36 1953]

then it cannot be decomposed until x is about 2/3, far far above the 0.0007 you have.

If you explore the form

[0.0007, -36;-36 x]

then it cannot be decomposed until x is about 1851400, which is about 950 times what you have, getting close to 1/2 the square of what you have.

So, your case is not borderline: your matrix appears to be firmly not positive definite.

Answer 2

[John BG]

Josh

your function E doesn't seem to work correctly for some input values, for instance

   E([5 -5])
  ans =
       0

but

E([1 1])
 =
   NaN
 E([1 1 1 1 1])
 =
   NaN

   mean([1 1 1 1 1])
   =
       1

if instead of E() mean is used, and taking into account that MATLAB definition of cov includes 1/(N-1) here it's 1/12 then replacing

cov = E((x-E(x)).*(z-E(z)));

with

cov = mean((x-mean(x)).*(z-mean(z)))*1/12

then the negative values that generate a negative eigenvalue is 2 orders of magnitude smaller

sigma =
   1.0e+03 *
   0.000000710897436  -0.000004496636588
  -0.000004496636588   1.952922555278590

then correcting (it error can be neglected) values <0 to 0:

sigma_ =
   1.0e+03 *
   0.000000710897436                   0
                   0   1.952922555278590

Y = mvncdf([0,0],X,mu,sigma_)

Y =
     6.151101794521690e-25

or perhaps you would like to consider abridging your code to the following:

clear all;

remove cov as variable otherwise MATLAB function cov doesn't work

x = [.125,.125,.125,.125,.125,.125,.125,.195,.195,.11,.135,.135,.145];
z = [466.087,480,483.4783,486.9565,466.087,424.3478,...
463.3043,468.2171,390.6977,335.7143,417.094,478.6325,450.4202];

mu=mean([x;z],2)
sigma=cov(x,z)
X = [1,1];
Y = mvncdf([0,0],X,mu',sigma)

Y =
     6.151102493733966e-25

it returns almost the same result

if you find this answer useful would you please be so kind to mark my answer as Accepted Answer?

To any other reader, please if you find this answer of any help solving your question,

please click on the thumbs-up vote link,

thanks in advance

John BG

<mailto:[e-mail] [e-mail]>