Friday, December 31, 2010

Generación de variables normales con una matriz de covarianzas

Este es un problema que tuve hace un tiempo, que lo solucione de otra forma (no me acuerdo como). Hace poco se me apareció esta solución, y me pareció interesante.
Nosotros tenemos una matriz de NxN de covarianzas (C). Queremos generar m samples aleatorios de N normales que respeten esa matriz de covarianzas.
Como la matriz C es real y simétrica, podemos encontrar una eigen decomposition:
C = RDR^tr
donde R es una matriz orthogonal formada por los autovectores, y D una matriz diagonal con los autovalores en la diagonal.
Por otro lado, si tenemos m realizaciones de n variables normales aleatorias independientes (con media cero) en una matriz X de mxn, la estimación de la matriz de Covarianzas sería:
C = 1/mX^trX.
Si las n variables son independientes, entonces, 1/mX^trX. va a ser una matriz diagonal, y el elemento {i,i} de la diagonal va a corresponder con la varianza de la variable i-esima.
Ahora, si generamos m realizaciones de n variables aleatorias, y cada variable aleatoria tiene media cero, y varianza igual al i-esimo autovalor, tenemos que
D = 1/mX^trX
Por lo que, si definimos a
U = XR^tr
tenemos una matriz de n variables con m realizaciones, y con matriz de covarianzas igual a
= 1/mU^trU
= 1/mRX^trXR^tr
= R( 1/mX^trX)R^tr
= RDR^tr
C = RDR^tr.

Que era lo que estábamos (bue, estaba buscando yo, no se si todos estábamos buscando eso) buscando.

Este metodo es super costoso si la matriz de covarianzas cambia siempre -por lo costoso de calcular los autovalores-, pero realmente sirve para entender por que PCA funciona. En realidad, habría que tener en cuenta los residuos, espero en algún momento del próximo año hacer un post sobre eso
Feliz año