Generación de muestras normales multivariantes: ¿por qué Cholesky?

Question

¡Hola a todos y feliz año nuevo! Que todas vuestras esperanzas y aspiraciones se hagan realidad y que las fuerzas del mal se vean confundidas y desorientadas de camino a vuestra casa.

Con eso fuera del camino...

Estoy intentando escribir un código informático que obtenga un vector $\mu \in R^n$ y la matriz $\Sigma \in \mathbb R^{n \times n}$ y genera muestras aleatorias de la distribución normal multivariante con media $\mu$ y la covarianza $\Sigma$ .

El problema : Sólo se me permite utilizar el programa para tomar muestras del solo distribución normal variable con media $0$ y la varianza $1$ : $N(0, 1)$ .

La solución propuesta : Definir un vector de ceros (inicialmente) $v \in \mathbb R^n$ , ahora para todos $i$ de $1$ à $n$ , extraer de una dist. normal de una sola variable: $v_i \overset{}{\sim} N(0, 1)$ .

Ahora haz una descomposición de Cholesky en $\Sigma$ : $\Sigma = LL^T$ .

Ahora finalmente el vector aleatorio que queremos que se distribuya de la gaussiana multivariante es $Lv + \mu$ .

Mi pregunta es ¿Por qué? No entiendo la intuición, si fuera una distribución unidimensional $N(\mu, \sigma^2)$ entonces entiendo por qué $\sigma ^2 v + \mu$ es una buena idea, ¿por qué Cholesky? ¿No querríamos $\Sigma v + \mu$ ?

Answer 1

Tras el comentario de Rahul has entendido que en cualquier parametrización $x=Av+μ$ necesitarás que $$Σ=\Bbb E(x-μ)(x-μ)^T=A·\Bbb E(vv^T)·A^T=AA^T.$$ Hay infinitas posibilidades para elegir $A$ con cualquier matriz ortogonal $Q$ también $\tilde A=AQ$ satisface esa condición.

Incluso se podría elegir la raíz cuadrada de $Σ$ (que existe y es única entre las matrices s.p.d.).

La ventaja de utilizar la factorización de Cholesky es que se dispone de un algoritmo barato y sencillo para calcularla.

Answer 2

Si todas las variables de la gaussiana multivariante fueran independientes, no tendríamos más problema que utilizar la fórmula $X_i =\sigma_i \nu+\mu_i$ . Como están correlacionados, tenemos (por ejemplo, el caso bivariante), $X_1 = \sigma_1\nu_1+\mu_1$ y $X_2 = \sigma_2[\rho\nu_1+\sqrt{1-\rho_{12}^2}\nu_2]+\mu_2$ y puede extenderse aún más a N. Nota: $$\sum = \begin{bmatrix}\sigma_1^2 &\rho_{12} \sigma_1\sigma_2 &\rho_{13} \sigma_2\sigma_3 & \dots \\\\\rho_{12} \sigma_1\sigma_2 &\sigma_2^2 &\rho_{23} \sigma_2\sigma_3 & \dots\end{bmatrix}$$ Al descomponer a través de Cholesky $\sum=LL^T$ podemos conseguir nuestro $X = L\nu+\mu$ sin cálculos manuales que, por otra parte, son bastante tediosos para el orden superior.