33
$\begingroup$

Bug introduced in 8.0 and fixed in 11.1

I found a strange behavior regarding the CDF of the bivariate Normal distribution

CDF[MultinormalDistribution[{0,0},({{1,37/40},{37/40,1}})],{0,0.2}]

gives

0.446357

On the other hand, the direct way

NIntegrate[PDF[MultinormalDistribution[{0,0},({{1,37/40},{37/40,1}})],{x, y}],{x,-\[Infinity],0},{y,-\[Infinity],0.2}]

gives

0.470073

which is right, confirmed by R software and other programs.

What the hell is going on here?

Improve this question
$\endgroup$
23
  • 4
    $\begingroup$ Please do not add the bugs tag when you post a new question. It is a special tag that is meant to be added by someone else than the original asker. $\endgroup$
    – Szabolcs
    Oct 29, 2016 at 12:53
  • 3
    $\begingroup$ If you look at this plot, those cusps shouldn't be there. They should be "round". ContourPlot[ CDF[d, {x, y}], {x, -3, 3}, {y, -3, 3}, PlotRange -> All, MaxRecursion -> 3 ]. Does look like a bug. Please do report this to Wolfram Support. $\endgroup$
    – Szabolcs
    Oct 29, 2016 at 13:17
  • 5
    $\begingroup$ @fwgb Hamming said the purpose of computing is insight not numbers. You have to learn to mistrust computers too. $\endgroup$
    – bobbym
    Oct 29, 2016 at 13:30
  • 6
    $\begingroup$ People who are voting to close as "simple mistake": what exactly is the mistake here? I don't see why this is the OP's mistake. It looks like a bug to me. A very disturbing bug that could invalidate results, as the OP said. I don't see good evidence that this is really a precision loss issue. Maybe it is, but I am not convinced. $\endgroup$
    – Szabolcs
    Oct 29, 2016 at 15:07
  • 4
    $\begingroup$ @fwgb Thank you for the reporting the issue. It has now been addressed in the sources and the fix will become available in the next release. The bug affects machine precision computation for $\rho^2 > 0.95^2$. Please use computation at $MachinePrecision as a temporary work-around. $\endgroup$
    – Sasha
    Nov 17, 2016 at 14:59

2 Answers 2

Reset to default
30
$\begingroup$

I almost believe the precision argument. But not quite.

dist = MultinormalDistribution[{0, 0}, ({{1, 37/40}, {37/40, 1}})];
CDF[dist, {0, 0.2`3}]
Precision[0.2]
$MachinePrecision
CDF[dist, {0, 0.2}]
(*  0.47  *)
(*  MachinePrecision  *)
(*  15.9546  *)
(*  0.446357  *)

So with only three digits of initial and intermediate precision, we get the right answer, but with nearly 16 (initially), we do not. Even assuming less than one digit of precision in the input, the correct output cannot be produced by the MachinePrecision computation. (The function is monotonic over the interval used.)

NMinimize[{CDF[dist, {0, y}], 0.0 <= y <= 0.4}, y]
NMaximize[{CDF[dist, {0, y}], 0.0 <= y <= 0.4}, y]
(*  {0.437977, {y -> 0.}}  *)
(*  {0.465287, {y -> 0.4}}  *)

The discrepancy between the correct CDF (blue) and the MachinePrecision CDF (yellow) can be quite large.

Plot[{
    NIntegrate[ PDF[ MultinormalDistribution[
      {0, 0}, ({{1, 37/40}, {37/40, 1}})], {x, y}], 
      {x, -∞, 0}, {y, -∞, u}],
    CDF[dist, {0, u}]}, 
  {u, 0, 1}]

Mathematica graphics

(I've weakly checked that the large discrepancy is not a result of numerical integration. Taking $m$ to be the off-diagonal element of the covariance matrix and assuming $0 < m < 1$, either the $x$ integral or the $y$ integral can be performed by Integrate, giving $\frac{1}{2 \sqrt{2 \pi}} \mathrm{e}^{-\frac{y^2}{2}} \text{erfc}\left(\frac{m y}{\sqrt{2-2 m^2}}\right)$ and $\frac{1}{2 \sqrt{2 \pi }} \mathrm{e}^{-\frac{x^2}{2}} \left(\text{erf}\left(\frac{u-m x}{\sqrt{2-2 m^2}}\right)+1\right)$, respectively. Replacing the double numerical integral with the single numerical integral of either of these does not visibly change the graph.)

Also, Karsten 7. is correct. This discrepancy suddenly turns on for a critical value of the off-diagonal covariance elements near 0.925.

Plot[
  CDF[ MultinormalDistribution[
      {0, 0}, ({{1, SetPrecision[x, 15]}, {SetPrecision[x, 15], 1}})], 
    {0, 0.2`15}] - 
  CDF[MultinormalDistribution[{0, 0}, ({{1, x}, {x, 1}})], 
    {0, 0.2}], 
  {x, 0.8, 1}, PlotRange -> All]

Mathematica graphics

This last is very strong evidence of a method switch introducing error, not precision loss.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thanks a lot. I was not aware that a different method is used for diagonal values larger than 0.925. Now the gap makes sense and only this method is corrupted. $\endgroup$
    – JHT
    Oct 30, 2016 at 9:33
11
$\begingroup$ 
dist = MultinormalDistribution[{0, 0}, ({{1, 37/40}, {37/40, 1}})];

MachinePrecision: "Machine-precision numbers (often called simply 'machine numbers') always contain a fixed number of digits and maintain no information about precision. ... Machine-precision computations are typically performed using native floating-point unit and low-level numeric library operations that are typically very fast (particularly so in matrix arithmetic), but provide no tracking of precision loss that may occur due to numerical round-off and other factors during a computation. As a result, machine arithmetic gives fast but numerically unvalidated results that may differ substantially from correct values."

cdf1 = CDF[dist, {0, .2}]

(*  0.446357  *)

Arbitrary-Precision Numbers "When you do calculations with arbitrary-precision numbers, the Wolfram Language keeps track of precision at all points. In general, the Wolfram Language tries to give you results which have the highest possible precision, given the precision of the input you provided."

Changing 0.2 to an arbitrary-precision number

cdf2 = CDF[dist, {0, .2`15}]

(*  0.47007303261424  *)

Similarly, using WorkingPrecision when plotting

ContourPlot[
 CDF[dist, {x, y}], {x, -3, 3}, {y, -3, 3},
 PlotRange -> All,
 MaxRecursion -> 3,
 WorkingPrecision -> 15]

enter image description here

Improve this answer
$\endgroup$
6
  • $\begingroup$ Good explanation, but changing 37/40 to 37./40. makes the above plot wrong again. How to aviod that? $\endgroup$
    – JHT
    Oct 29, 2016 at 14:58
  • $\begingroup$ @fwgb Make sure you do not introduce any machine precision numbers. 37. is machine precision. MachinePrecision is contagious and will turn everything else it touches into machine precision. SetPrecision can convert to arbitrary precision. Rationalize can convert to exact. $\endgroup$
    – Szabolcs
    Oct 29, 2016 at 15:05
  • 1
    $\begingroup$ I don't think that is the full story here. There is something else going wrong. I suspect that the implementation switches to a buggy implementation, when one parameter has MachinePrecision. $\endgroup$
    – Karsten7
    Oct 29, 2016 at 17:49
  • 1
    $\begingroup$ @Karsten7. Yes you're right. The funny thing is that for <= 36/40 everything is fine. For >= 37/40 the problem occurs. $\endgroup$
    – JHT
    Oct 29, 2016 at 18:36
  • 7
    $\begingroup$ @fwgb Glancing over the code, it looks like Statistics`MultinormalDistributionsDump`GenzDreznerWeselowskyBvnlCDF gets called, which contains an If[ar < 0.925, ...]. Test: Plot[CDF[MultinormalDistribution[{0, 0}, ({{1, SetPrecision[x, 15]}, {SetPrecision[x, 15], 1}})], {0, 0.2`15}] - CDF[MultinormalDistribution[{0, 0}, ({{1, x}, {x, 1}})], {0, 0.2}], {x, 0.8, 1}, PlotRange -> All] $\endgroup$
    – Karsten7
    Oct 29, 2016 at 20:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.