Serious problem with CDF[MultinormalDistribution...
Myeong Ae Kang
I have a problem, serious and urgent.
I got wrong results of CDF of MultinormalDistribution.
I'm calculating (1.0 -CDF[MultinormalDistribution[{0, 0,0}, {{1, 0.559887,
0.289988}, {0.559887, 1,0.880157}, {0.289988, 0.880157, 1}}],
{4.5,4.5,6.9958}]).
The results is smaller than (1.0 -CDF[MultinormalDistribution[{0,0}, {{1,
0.559887}, {0.559887, 1}}], {4.5,4.5}]). It should be larger.
Actually, I have to calculate more than 6 to 12 random variables (12 x 12 sigma
matrix.).
Mathematica 6.0 and 7.0.1 gave me the same wrong result.
I compared this with the result of another system.
That system gave the expected value.
To fix this problem, what should I do?
Regards,
Myeong Ae Kang
Darren Glosemeyer
This is addressed in version 8. Here are the results:
In[1]:= threeD = (1.0 -
CDF[MultinormalDistribution[{0, 0,
0}, {{1, 0.559887, 0.289988}, {0.559887, 1,
0.880157}, {0.289988, 0.880157, 1}}], {4.5, 4.5, 6.9958}])
-6
Out[1]= 6.75155 10
In[2]:= twoD = (1.0 -
CDF[MultinormalDistribution[{0,
0}, {{1, 0.559887}, {0.559887, 1}}], {4.5, 4.5}])
-6
Out[2]= 6.75155 10
In[3]:= threeD > twoD
Out[3]= True
The difference is at the level of machine precision numerical error:
In[4]:= threeD - twoD
-16
Out[4]= 1.11022 10
Note that if we artificially increase the precision as in the example
below, we get a difference at the same order of magnitude, so the
precision in the difference above is as good as could be expected (and
in fact probably better than should be expected) at machine precision.
In[5]:= CDF[MultinormalDistribution[{0, 0, 0},
Rationalize[{{1, 0.559887, 0.289988}, {0.559887, 1,
0.880157}, {0.289988, 0.880157, 1}}, 0]],
SetPrecision[{4.5, 4.5, 6.9958}, 50]]
Out[5]= 0.999993248446123977600817893853727681422304494350
In[6]:= CDF[MultinormalDistribution[{0, 0},
Rationalize[{{1, 0.559887}, {0.559887, 1}}, 0]],
SetPrecision[{4.5, 4.5}, 50]]
Out[6]= 0.999993248446124136306526964076086153533229885505
In[7]:= % - %%
-16
Out[7]= 1.5870570907022235847211092539116 10
Darren Glosemeyer
Wolfram Research
Myeong Ae Kang
values, the results are much worse even than this.
In[1326]:= Needs["MultivariateStatistics`"]
(1.0-CDF[MultinormalDistribution[{0,0,0},Rationalize[{{1,0.559887,0.289988},
{0.559887,1,0.880157},{0.289988,0.880157,1}}]],SetPrecision
[{4.5,4.5,6.9958},50]])
(1.0-CDF[MultinormalDistribution[{0,0},{{1,0.559887},{0.559887,1}}],{4.5,4.5}])
Out[1327]= 5.00832*10^-6
Out[1328]= 6.75155*10^-6
And, similarly to this, with version 6.0, I did the calculation of CDF with 6
RV's in different orders of RV's. I just changed the orders of the sigma matrix
and range vector with zero mean. But, the results(which are suppossed to equal)
are quite different like this following:
In[1432]:= ClearAll;
Print["R,A,W ======================="]
TT={{1,0.576363,0.300013,0.888477,0.106961,0.215919},
{0.576363,1,0.888477,0.300013,0.215919,0.106961},{0.300013,0.888477,1,-
9.72249*10^-6,0.466809,-0.191847},{0.888477,0.300013,-9.72249*10^-6,1,-
0.191847,0.466809},{0.106961,0.215919,0.466809,-0.191847,1,-0.862483},
{0.215919,0.106961,-0.191847,0.466809,-0.862483,1}};
TTbnd={7.14659,7.14659,7.83599,7.83599,5.55663,5.55663};
TTp={{0,0,0,0,0,1},{0,0,0,0,1,0},{0,0,0,1,0,0},{0,0,1,0,0,0},{0,1,0,0,0,0},
{1,0,0,0,0,0}}.TT.{{0,0,0,0,0,1},{0,0,0,0,1,0},{0,0,0,1,0,0},{0,0,1,0,0,0},
{0,1,0,0,0,0},{1,0,0,0,0,0}}
Tbnd=TTbnd.Transpose[{{0,0,0,0,0,1},{0,0,0,0,1,0},{0,0,0,1,0,0},{0,0,1,0,0,0},
{0,1,0,0,0,0},{1,0,0,0,0,0}}]
1.0-CDF[MultinormalDistribution[{0,0,0,0,0,0},TTp],Tbnd]
Print["A,R,W ======================="]
TTp2={{0,0,1,0,0,0},{0,0,0,1,0,0},{0,0,0,0,1,0},{0,0,0,0,0,1},{1,0,0,0,0,0},
{0,1,0,0,0,0}}.TT.Transpose[{{0,0,1,0,0,0},{0,0,0,1,0,0},{0,0,0,0,1,0},
{0,0,0,0,0,1},{1,0,0,0,0,0},{0,1,0,0,0,0}}]
Tbnd2=TTbnd.Transpose[{{0,0,1,0,0,0},{0,0,0,1,0,0},{0,0,0,0,1,0},{0,0,0,0,0,1},
{1,0,0,0,0,0},{0,1,0,0,0,0}}]
1.0-CDF[MultinormalDistribution[{0,0,0,0,0,0},TTp2],Tbnd2]
Print["A,W,R ======================="]
TTp3={{0,0,1,0,0,0},{0,0,0,1,0,0},{1,0,0,0,0,0},{0,1,0,0,0,0},{0,0,0,0,1,0},
{0,0,0,0,0,1}}.TT.Transpose[{{0,0,1,0,0,0},{0,0,0,1,0,0},{1,0,0,0,0,0},
{0,1,0,0,0,0},{0,0,0,0,1,0},{0,0,0,0,0,1}}]
Tbnd3=TTbnd.Transpose[{{0,0,1,0,0,0},{0,0,0,1,0,0},{1,0,0,0,0,0},{0,1,0,0,0,0},
{0,0,0,0,1,0},{0,0,0,0,0,1}}]
1.0-CDF[MultinormalDistribution[{0,0,0,0,0,0},TTp3],Tbnd3]
During evaluation of In[1432]:= R,A,W =======================
Out[1436]= {{1.,-0.862483,0.466809,-0.191847,0.106961,0.215919},{-0.862483,1.,-
0.191847,0.466809,0.215919,0.106961},{0.466809,-0.191847,1.,-9.72249*10^-
6,0.300013,0.888477},{-0.191847,0.466809,-9.72249*10^-6,1.,0.888477,0.300013},
{0.106961,0.215919,0.300013,0.888477,1.,0.576363},
{0.215919,0.106961,0.888477,0.300013,0.576363,1.}}
Out[1437]= {5.55663,5.55663,7.83599,7.83599,7.14659,7.14659}
Out[1438]= 1.37516*10^-8
During evaluation of In[1432]:= A,R,W =======================
Out[1440]= {{1.,-9.72249*10^-6,0.466809,-0.191847,0.300013,0.888477},{-
9.72249*10^-6,1.,-0.191847,0.466809,0.888477,0.300013},{0.466809,-0.191847,1.,-
0.862483,0.106961,0.215919},{-0.191847,0.466809,-0.862483,1.,0.215919,0.106961},
{0.300013,0.888477,0.106961,0.215919,1.,0.576363},
{0.888477,0.300013,0.215919,0.106961,0.576363,1.}}
Out[1441]= {7.83599,7.83599,5.55663,5.55663,7.14659,7.14659}
Out[1442]= 6.50208*10^-9
During evaluation of In[1432]:= A,W,R =======================
Out[1444]= {{1.,-9.72249*10^-6,0.300013,0.888477,0.466809,-0.191847},{-
9.72249*10^-6,1.,0.888477,0.300013,-0.191847,0.466809},
{0.300013,0.888477,1.,0.576363,0.106961,0.215919},
{0.888477,0.300013,0.576363,1.,0.215919,0.106961},{0.466809,-
0.191847,0.106961,0.215919,1.,-0.862483},{-0.191847,0.466809,0.215919,0.106961,-
0.862483,1.}}
Out[1445]= {7.83599,7.83599,7.14659,7.14659,5.55663,5.55663}
Out[1446]= 1.17684*10^-14
---------------
Now, I have two questions.
1) Does mathematica version 8.0 give three similar results of the second part?
I don't have version 8.0. Before upgrading it (if possible), I wonder that.
2) Is there any other way to get the proper results with version 6.0 or 7.01?
Thank you.
Regards,
Myeong Ae Kang
=================================