Cómo utilizar la descomposición de Cholesky, o una alternativa, para la simulación de datos correlacionados

Question

Utilizo la descomposición de Cholesky para simular variables aleatorias correlacionadas dada una matriz de correlaciones. El caso es que el resultado nunca reproduce la estructura de correlación tal y como se da. He aquí un pequeño ejemplo en Python para ilustrar la situación.

import numpy as np    

n_obs = 10000
means = [1, 2, 3]
sds = [1, 2, 3] # standard deviations 

# generating random independent variables 
observations = np.vstack([np.random.normal(loc=mean, scale=sd, size=n_obs)
                   for mean, sd in zip(means, sds)])  # observations, a row per variable

cor_matrix = np.array([[1.0, 0.6, 0.9],
                       [0.6, 1.0, 0.5],
                       [0.9, 0.5, 1.0]])

L = np.linalg.cholesky(cor_matrix)

print(np.corrcoef(L.dot(observations))) 

Esto imprime:

[[ 1.          0.34450587  0.57515737]
 [ 0.34450587  1.          0.1488504 ]
 [ 0.57515737  0.1488504   1.        ]]

Como puede ver, la matriz de correlaciones estimada post hoc difiere drásticamente de la anterior. Hay algún error en mi código o existe alguna alternativa al uso de la descomposición Cholesky?

Editar

Le pido perdón por este lío. No creía que hubiera un error en el código y/o en la forma de aplicar la descomposición de Cholesky debido a algún malentendido del material que había estudiado antes. De hecho estaba seguro de que el método en sí no pretendía ser preciso y me había parecido bien hasta la situación que me hizo publicar esta pregunta. Gracias por señalarme el error de concepto que tenía. He editado el título para reflejar mejor la situación real propuesta por @Silverfish.

Answer 1

Como ya han demostrado otros: cholesky funciona. Aquí un trozo de código que es muy corto y muy cerca de pseudocódigo: un codepiece en MatMate:

Co = {{1.0, 0.6, 0.9},  _
      {0.6, 1.0, 0.5},  _
      {0.9, 0.5, 1.0}}           // make correlation matrix

chol = cholesky(co)              // do cholesky-decomposition           
data = chol * unkorrzl(randomn(3,100,0,1))  
                                 // dot-multiply cholesky with random-
                                 // vectors with mean=0, sdev=1  
                                 //(refined by a "decorrelation" 
                                 //to remove spurious/random correlations)   

chk = data *' /100               // check the correlation of the data
list chk

1.0000  0.6000  0.9000
0.6000  1.0000  0.5000
0.9000  0.5000  1.0000