Necesito aplicar scipy.stats.multivariate_normal.cdf(), que calcula la 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$$ ,
donde $\mu$ es la media y $\Sigma^{-1}$ es la matriz inversa de la matriz de covarianza $\Sigma$ .
Para que me entienda, estoy intentando realizar el mismo cálculo en otro programa (Maple). Sea
\Sigma = \begin{pmatrix} a_{11} & a_{12} & a_{13}\\ a_{21} & a_{22} & a_{23}\\ a_{31} & a_{32} & a_{33} \end{pmatrix}
y a efectos de este cálculo dejemos que $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$ .
Si expandimos el término de la exponencial, obtenemos
$$-\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)}$$
y así la integral que estoy tratando de evaluar en Maple es:
$$\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$$ .
El resultado es $4.164024864\times10^{-242}$ .
Aplicando scipy.stats.multivariate_normal.cdf() con $[x, y, z] = [1.37824, -21.58961, 18.48617]$ , $mean=None$ y $cov=np.array([[1, 0.87055, 0.710804075], [0.87055, 1, 0.8165], [0.710804075, 0.8165, 1]])$ .
El resultado enPython es $1.124512788731174 \times 10^{-103}$ .
La diferencia entre los dos resultados es significativa. Parece que estoy haciendo algo mal. Agradeceré si alguien puede señalar dónde está mi error. Gracias de antemano.
