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!
