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

조회 수: 6 (최근 30일)
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
Walter Roberson 2017년 3월 9일
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.

추가 답변 (1개)

John BG
John BG 2017년 3월 9일
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

카테고리

Help CenterFile Exchange에서 Detection에 대해 자세히 알아보기

제품

Community Treasure Hunt

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

Start Hunting!

Translated by