30
$\begingroup$

Bug introduced in 8.0 or earlier 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?

$\endgroup$
  • 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 '16 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 '16 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 '16 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 '16 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 '16 at 14:59
29
$\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.

$\endgroup$
  • $\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$ – fwgb Oct 30 '16 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

$\endgroup$
  • $\begingroup$ Good explanation, but changing 37/40 to 37./40. makes the above plot wrong again. How to aviod that? $\endgroup$ – fwgb Oct 29 '16 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 '16 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$ – Karsten 7. Oct 29 '16 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$ – fwgb Oct 29 '16 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$ – Karsten 7. Oct 29 '16 at 20:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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