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

jg614

### jg614 (view profile)

님이 질문을 제출함. 9 Mar 2017
최근 활동 John BG

### John BG (view profile)

님이 답변함. 9 Mar 2017
Walter Roberson

### Walter Roberson (view profile)

님이 답변을 채택함.
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.

로그인 to comment.

## 답변 수: 2

### Walter Roberson (view profile)

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.

로그인 to comment.

### John BG (view profile)

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