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

Utilizo la descomposición de Cholesky para simular variables aleatorias correlacionadas dada una matriz de correlación. La cuestión es que el resultado nunca reproduce la estructura de correlación tal como se da. Aquí hay 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 correlación estimada post-hoc difiere drásticamente de la anterior. ¿Hay algún error en mi código o hay alguna alternativa para usar la descomposición de Cholesky?

Editar

Perdón por este desastre. No pensé que hubiera un error en el código y / o en la forma en que se aplicó la descomposición de Cholesky debido a algunos malentendidos del material que había estudiado antes. De hecho, estaba seguro de que el método en sí no estaba destinado a ser preciso y había estado de acuerdo con eso hasta la situación que me hizo publicar esta pregunta. Gracias por señalar el error que tuve. He editado el título para reflejar mejor la situación real propuesta por @Silverfish.

El enfoque basado en la descomposición de Cholesky debería funcionar, se describe aquí y se muestra en la respuesta de Mark L. Stone publicado casi al mismo tiempo que esta respuesta.

$N(\vec\mu, \Sigma)$

$Y$$X$$\Phi$$\Sigma$$\Lambda$$\Sigma$$\Phi$

Ejemplo en R(lo siento, no estoy usando el mismo software que usaste en la pregunta):

n <- 10000
corM <- rbind(c(1.0, 0.6, 0.9), c(0.6, 1.0, 0.5), c(0.9, 0.5, 1.0))
set.seed(123)
SigmaEV <- eigen(corM)
eps <- rnorm(n * ncol(SigmaEV$vectors))
Meps <- matrix(eps, ncol = n, byrow = TRUE)    
Meps <- SigmaEV$vectors %*% diag(sqrt(SigmaEV$values)) %*% Meps
Meps <- t(Meps)
# target correlation matrix
corM
#      [,1] [,2] [,3]
# [1,]  1.0  0.6  0.9
# [2,]  0.6  1.0  0.5
# [3,]  0.9  0.5  1.0
# correlation matrix for simulated data
cor(Meps)
#           [,1]      [,2]      [,3]
# [1,] 1.0000000 0.6002078 0.8994329
# [2,] 0.6002078 1.0000000 0.5006346
# [3,] 0.8994329 0.5006346 1.0000000

También te puede interesar esta publicación y esta publicación .

Es probable que la gente encuentre su error mucho más rápido si explica lo que hizo con palabras y álgebra en lugar de código (o al menos escribirlo con pseudocódigo).

Parece que estás haciendo el equivalente de esto (aunque posiblemente transpuesto):

1. $n\times k$$Z$

2. $\sigma_i$$\mu_i$

3. $Y=LX$

$L$

Lo que debes hacer es esto:

1. $n\times k$$Z$

2. $X=LZ$

3. $\sigma_i$$\mu_i$

Hay muchas explicaciones de este algoritmo en el sitio. p.ej

¿Cómo generar números aleatorios correlacionados (medias, variaciones y grado de correlación dados)?

¿Puedo usar el método Cholesky para generar variables aleatorias correlacionadas con la media dada?

Éste lo discute directamente en términos de la matriz de covarianza deseada, y también proporciona un algoritmo para obtener una covarianza de muestra deseada :

Generando datos con una matriz de covarianza de muestra dada

No hay nada malo con la factorización de Cholesky. Hay un error en su código. Ver edición a continuación.

Aquí está el código y los resultados de MATLAB, primero para n_obs = 10000 como lo ha hecho, luego para n_obs = 1e8. Para simplificar, ya que no afecta los resultados, no me molesto con los medios, es decir, los hago ceros. Tenga en cuenta que el chol de MATLAB produce un factor de Cholesky triangular superior R de la matriz M tal que R '* R = M. numpy.linalg.cholesky produce un factor de Cholesky triangular inferior, por lo que se necesita un ajuste en función de mi código; pero creo que su código está bien en ese sentido.

   >> correlation_matrix = [1.0, 0.6, 0.9; 0.6, 1.0, 0.5;0.9, 0.5, 1.0];
   >> SD = diag([1 2 3]);
   >> covariance_matrix = SD*correlation_matrix*SD
   covariance_matrix =
      1.000000000000000   1.200000000000000   2.700000000000000
      1.200000000000000   4.000000000000000   3.000000000000000
      2.700000000000000   3.000000000000000   9.000000000000000
   >> n_obs = 10000;
   >> Random_sample = randn(n_obs,3)*chol(covariance_matrix);
   >> disp(corr(Random_sample))
      1.000000000000000   0.599105015695768   0.898395949647890
      0.599105015695768   1.000000000000000   0.495147514173305
      0.898395949647890   0.495147514173305   1.000000000000000
   >> n_obs = 1e8;
   >> Random_sample = randn(n_obs,3)*chol(covariance_matrix);
   >> disp(corr(Random_sample))
      1.000000000000000   0.600101477583914   0.899986072541418
      0.600101477583914   1.000000000000000   0.500112824962378
      0.899986072541418   0.500112824962378   1.000000000000000

Editar: encontré tu error. Usted aplicó incorrectamente la desviación estándar. Esto es el equivalente de lo que hiciste, lo cual está mal.

   >> n_obs = 10000;
   >> Random_sample = randn(n_obs,3)*SD*chol(correlation_matrix);
   >> disp(corr(Random_sample))
      1.000000000000000   0.336292731308138   0.562331469857830
      0.336292731308138   1.000000000000000   0.131270077244625
      0.562331469857830   0.131270077244625   1.000000000000000
   >> n_obs=1e8;
   >> Random_sample = randn(n_obs,3)*SD*chol(correlation_matrix);
   >> disp(corr(Random_sample))
      1.000000000000000   0.351254525742470   0.568291702131030
      0.351254525742470   1.000000000000000   0.140443281045496
      0.568291702131030   0.140443281045496   1.000000000000000

El CV no se trata de código, pero me intrigó ver cómo se vería esto después de todas las buenas respuestas, y específicamente la contribución de @Mark L. Stone. La respuesta real a la pregunta se proporciona en su publicación (acredite su publicación en caso de duda). Estoy moviendo esta información adjunta aquí para facilitar la recuperación de esta publicación en el futuro. Sin restarle importancia a ninguna de las otras respuestas excelentes, después de la respuesta de Mark, esto cierra el problema corrigiendo la publicación en el OP.

Fuente

En Python:

import numpy as np

no_obs = 1000             # Number of observations per column
means = [1, 2, 3]         # Mean values of each column
no_cols = 3               # Number of columns

sds = [1, 2, 3]           # SD of each column
sd = np.diag(sds)         # SD in a diagonal matrix for later operations

observations = np.random.normal(0, 1, (no_cols, no_obs)) # Rd draws N(0,1) in [3 x 1,000]

cor_matrix = np.array([[1.0, 0.6, 0.9],
                       [0.6, 1.0, 0.5],
                       [0.9, 0.5, 1.0]])          # The correlation matrix [3 x 3]

cov_matrix = np.dot(sd, np.dot(cor_matrix, sd))   # The covariance matrix

Chol = np.linalg.cholesky(cov_matrix)             # Cholesky decomposition

array([[ 1.        ,  0.        ,  0.        ],
       [ 1.2       ,  1.6       ,  0.        ],
       [ 2.7       , -0.15      ,  1.29903811]])

sam_eq_mean = Chol .dot(observations)             # Generating random MVN (0, cov_matrix)

s = sam_eq_mean.transpose() + means               # Adding the means column wise
samples = s.transpose()                           # Transposing back

print(np.corrcoef(samples))                       # Checking correlation consistency.

[[ 1.          0.59167434  0.90182308]
 [ 0.59167434  1.          0.49279316]
 [ 0.90182308  0.49279316  1.        ]]

EN [R]:

no_obs = 1000             # Number of observations per column
means = 1:3               # Mean values of each column
no_cols = 3               # Number of columns

sds = 1:3                 # SD of each column
sd = diag(sds)         # SD in a diagonal matrix for later operations

observations = matrix(rnorm(no_cols * no_obs), nrow = no_cols) # Rd draws N(0,1)

cor_matrix = matrix(c(1.0, 0.6, 0.9,
                      0.6, 1.0, 0.5,
                      0.9, 0.5, 1.0), byrow = T, nrow = 3)     # cor matrix [3 x 3]

cov_matrix = sd %*% cor_matrix %*% sd                          # The covariance matrix

Chol = chol(cov_matrix)                                        # Cholesky decomposition

     [,1] [,2]      [,3]
[1,]    1  1.2  2.700000
[2,]    0  1.6 -0.150000
[3,]    0  0.0  1.299038

sam_eq_mean = t(observations) %*% Chol          # Generating random MVN (0, cov_matrix)

samples = t(sam_eq_mean) + means

cor(t(samples))

          [,1]      [,2]      [,3]
[1,] 1.0000000 0.6071067 0.8857339
[2,] 0.6071067 1.0000000 0.4655579
[3,] 0.8857339 0.4655579 1.0000000

colMeans(t(samples))
[1] 1.035056 2.099352 3.065797
apply(t(samples), 2, sd)
[1] 0.9543873 1.9788250 2.8903964

Como otros ya han demostrado: cholesky funciona. Aquí una pieza de código que es muy corta y muy cercana al pseudocódigo: una pieza de código 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