La pregunta que realmente le preocupa a los pares de la normal de variables. Vamos a llamarlos $x_1$$x_2$, con los medios de $\mu_i$, desviaciones estándar $\sigma_i$, y la correlación $\rho$. De donde sus conjunta pdf
$$\frac{1}{2 \pi \sqrt{1 - \rho^2} \sigma_1 \sigma_2}
e^{-\frac{1}{1-\rho^2} \left(\frac{(x_1 - \mu_1)^2}{2 \sigma_1^2} + \frac{(x_2 - \mu_2)^2}{2 \sigma_2^2} - \frac{\rho (x_1 - \mu_1)(x_2 - \mu_2)}{\sigma_1 \sigma_2}\right)} dx_1 dx_2\text{.}$$
Deje $f(x_1,x_2)$ ser el producto de este con el $y_i$ (como funciones de la $x_i$). El primer componente del gradiente de $\log(f)$ es
$$\frac{\partial \log(f)}{\partial x_1}
= \frac{1}{1 + e^{x_1}} + \frac{\rho(\mu_2 - x_2) \sigma_1 + (x_1 - \mu_1)\sigma_2}{(\rho^2-1)\sigma_1^2 \sigma_2},$$
con una expresión similar para el segundo componente (a través de la simetría alcanzado mediante el intercambio de los subíndices 1 y 2). Habrá un único máximo global, que se puede detectar mediante el ajuste de la pendiente a cero. Este par de ecuaciones no lineales no tiene solución de forma cerrada. Se encontró rápidamente por un par de Newton-Raphson iteraciones. Alternativamente, se puede linealizar las ecuaciones. En efecto, a través de segundo orden, el primer componente es igual a
$$\frac{1}{2} + x_1\left(\frac{-1}{4} + \frac{1}{(\rho^2-1)\sigma_1^2}\right) + \frac{-\rho x_2 \sigma_1 + \rho \mu_2 \sigma_1 - \mu_1 \sigma_2}{(\rho^2 -1)\sigma_1^2 \sigma_2}.$$
Esto le da un par de lineal de ecuaciones en $(x_1, x_2)$, lo que, por tanto, no tiene una forma cerrada de la solución, decir $\hat{x}_i(\mu_1, \mu_2, \sigma_1, \sigma_2, \rho)$, que obviamente son racionales de polinomios.
El Jacobiano en este momento crítico tiene 1,1 coeficiente de
$$\frac{e^\hat{x_1}\left(2 - (\rho^2-1)\sigma_1^2 + 2\cosh(\hat{x_1})\right)}{(1+e^\hat{x_1})^2(\rho^2-1)\sigma_1^2},$$
1,2 y 2,1 coeficientes de
$$\frac{\rho}{\sigma_1 \sigma_2(1 - \rho^2)},$$
y 2,2 coeficiente obtenido del coeficiente 1,1 por la simetría. Debido a que este es un punto crítico (al menos aproximadamente), se puede sustituir
$$e^\hat{x_1} = \frac{(\rho^2-1)\sigma_1^2 \sigma_2}{(\mu_2 - \hat{x_2})\rho \sigma_1 + (\hat{x_1} - \mu_1)\sigma_2} - 1$$
y el uso que también para calcular los $\cosh(\hat{x_1}) = \frac{e^\hat{x_1} - e^{-\hat{x_1}}}{2}$, con una manipulación similar para$e^\hat{x_2}$$\cosh(\hat{x_2})$. Esto permite la evaluación de la Arpillera (el determinante del Jacobiano) como una función racional de los parámetros.
El resto es de rutina: el de Hesse nos dice cómo aproximar la integral como una binormal integral (una saddlepoint aproximación). La respuesta es igual a $\frac{1}{2\pi}$ veces una función racional de los cinco parámetros: esa es su forma cerrada (por lo que vale la pena!).
