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

27 views (last 30 days)
jg614 on 9 Mar 2017
Answered: John BG on 9 Mar 2017
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.

Walter Roberson on 9 Mar 2017
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.
[ 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.

John BG on 9 Mar 2017
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
John BG

### Community Treasure Hunt

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

Start Hunting!

Translated by