2
$\begingroup$

I need to apply scipy.stats.multivariate_normal.cdf(), which computes the integral

$$\int \frac{1}{\sqrt{(2\pi)^{\frac{3}{2}}\det(\Sigma)}} \exp\left(-\frac{1}{2}(x-\mu)^T \Sigma^{-1}(x-\mu)\right) dx$$,

where $\mu$ is the mean and $\Sigma^{-1}$ is the inverse matrix of the covariance matrix $\Sigma$.

For my own understanding I am trying to run the same calculation in another program (Maple). Let

\Sigma = \begin{pmatrix} a_{11} & a_{12} & a_{13}\\ a_{21} & a_{22} & a_{23}\\ a_{31} & a_{32} & a_{33} \end{pmatrix}

and for the purpose of this computation let $a_{11} = a_{22} = a_{33} = 1; a_{12} = a_{21} = 0.87055, a_{13} = a_{31} = 0.710804075, a_{23} = a_{32} = 0.8165$.

If we expand the term in the exponential, we get

$$-\frac{1}{2}(x-\mu)^T \Sigma^{-1}(x-\mu) = -\frac{x^2}{2} - \frac{(y - 0.87055 x)^2}{2(1 - 0.87055^2)} - \frac{(z - 0.8165 y)^2}{2(1 - 0.8165^2)}$$

and so the integral I am trying to evaluate in Maple is:

$$\int\limits_{-\infty}^{1.37824} \int\limits_{-\infty}^{-21.58961} \int\limits_{-\infty}^{18.48617} \frac{\exp\left( -\frac{x^2}{2} - \frac{(y - 0.87055 x)^2}{2(1 - 0.87055^2)} - \frac{(z - 0.8165 y)^2}{2(1 - 0.8165^2)} \right)}{(2\pi)^{\frac{3}{2}}\sqrt{(1-0.87055^2)(1-0.8165^2)}} dz dy dx$$.

The result is $4.164024864\times10^{-242}$.

Applying scipy.stats.multivariate_normal.cdf() with $[x, y, z] = [1.37824, -21.58961, 18.48617]$, $mean=None$ and $cov=np.array([[1, 0.87055, 0.710804075], [0.87055, 1, 0.8165], [0.710804075, 0.8165, 1]])$.

The result inPython is $1.124512788731174 \times 10^{-103}$.

The difference between the two results is significant. Apparently I am doing something wrong. I will appreciate if anyone can point where my mistake is! Thank you in advance!

$\endgroup$
7
  • 5
    $\begingroup$ The difference is really, really small? $\endgroup$
    – Bob Jansen
    Nov 27 '21 at 22:46
  • 5
    $\begingroup$ You should compare the two probabilities in a region where they have sufficient mass. That upper limit of -21.58961 is too extreme in the negative. Try a selection of $[x,y,z]$ where the variables are between $-3$ and $+3\,.$ I would also try different values of the $a_{ij}\,.$ $\endgroup$
    – Kurt G.
    Nov 28 '21 at 6:37
  • $\begingroup$ @Bob Jansen, the difference is infinitesimal, indeed. But my expectation is that the two outputs should be close. However, it turns out that the powers of 10 are significant. $\endgroup$
    – Candidate
    Nov 28 '21 at 8:56
  • $\begingroup$ @Kurt G, I tried with the standard case (which can be even computed by hand), i.e. when $[x, y, z] = [0, 0, 0]$ and when the covariance matrix is \begin{pmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end{pmatrix}, i.e., $a_{ij}$ is $1$, when $i=j$ and $0$ otherwise, then Python and Maple yield $\frac{1}{8}$, which seems to be as expected. That only boosted the impression that I am doing something wrong. In any case, thank you both for your help! $\endgroup$
    – Candidate
    Nov 28 '21 at 8:56
  • 3
    $\begingroup$ Ok. Python and Maple are spot on for 0,0,0 and uncorrelated. The extrem example in your question above leads to zero for Python and Maple. Every numercial practitioner knows that numbers like that are to be treated as zero. You will even get answers with larger differences of orders 1E-15 if you do the calculation on different types of CPUs. Don't worry be happy. $\endgroup$
    – Kurt G.
    Nov 28 '21 at 9:01
2
$\begingroup$

Both values are practically zero and the difference between the two pieces of software is really, really small (almost infinitesimal compared to your inputs). Computers are not perfect calculators and some rounding error is to be expected. Try to be aware that they can happen and check whether the rounding error matters for your application.

$\endgroup$
1
  • $\begingroup$ Thank you @Bob Jansen! You are totally right, now I realise that the two outputs don't necessarily need to be the same (although they are practically zero, as the comments above imply). $\endgroup$
    – Candidate
    Nov 30 '21 at 20:54

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.