# Triple Gaussian integrals in Python

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!

• The difference is really, really small? Nov 27 '21 at 22:46
• 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}\,.$ Nov 28 '21 at 6:37
• @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. Nov 28 '21 at 8:56
• @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! Nov 28 '21 at 8:56
• 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. Nov 28 '21 at 9:01