24 votos

Cómo utilizar 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 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.

19voto

AdamSane Puntos 1825

La gente probablemente encontraría tu error mucho más rápido si explicaras lo que hiciste con palabras y álgebra en lugar de código (o al menos escribiéndolo usando pseudocódigo).

Usted parece estar haciendo el equivalente de esto (aunque posiblemente transpuesto):

  1. Generar un $n\times k$ matriz de normales estándar, $Z$

  2. multiplicar las columnas por $\sigma_i$ y añada $\mu_i$ para obtener normales no estándar

  3. calcular $Y=LX$ para obtener normales correlacionadas.

donde $L$ es el factor Cholesky izquierdo de su matriz de correlación.

Lo que debes hacer es lo siguiente:

  1. Generar un $n\times k$ matriz de normales estándar, $Z$

  2. calcular $X=LZ$ para obtener normales correlacionadas.

  3. multiplicar las columnas por $\sigma_i$ y añada $\mu_i$ para obtener normales no estándar

Hay muchas explicaciones de este algoritmo en la web, por ejemplo

¿Cómo generar números aleatorios correlacionados (dadas las medias, las varianzas y el grado de correlación)?

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

Este lo discute directamente en términos de la matriz de covarianza deseada, y también da un algoritmo para obtener un muestra covarianza:

Generar datos con una matriz de covarianza muestral dada

12voto

einverne Puntos 126

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

No obstante, a veces he generado extracciones de la distribución Normal multivariante $N(\vec\mu, \Sigma)$ como sigue:

$$ Y = Q X + \vec\mu \,, \quad \hbox{with}\quad Q=\Lambda^{1/2}\Phi \,, $$

donde $Y$ son los empates finales, $X$ son extracciones de la distribución normal estándar univariante, $\Phi$ es una matriz que contiene los vectores propios normalizados de la matriz objetivo $\Sigma$ et $\Lambda$ es una matriz diagonal que contiene los valores propios de $\Sigma$ dispuestos en el mismo orden que los vectores propios en las columnas de $\Phi$ .

Ejemplo en R (lo siento, no estoy utilizando el mismo software que utilizó 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 le puede interesar esta entrada y esta entrada .

11voto

Mark L. Stone Puntos 2037

La factorización Cholesky no tiene nada de malo. Hay un error en tu código. Ver editar a continuación.

Aquí está el código MATLAB y los resultados, primero para n_obs = 10000 como usted tiene, entonces para n_obs = 1e8. Por simplicidad, ya que no afecta a los resultados, no me molesto con las medias, es decir, las hago ceros. Ten en cuenta que 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 es necesario un ajuste frente a mi código; pero creo que tu 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

Edición: He encontrado tu error. Has aplicado incorrectamente la desviación típica. Esto es equivalente a lo que hiciste, que es incorrecto.

   >> 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

9voto

Antoni Parellada Puntos 2762

CV no se trata de código, pero me intrigaba ver cómo quedaría después de todas las buenas respuestas, y en concreto la contribución de @Mark L. Stone. La respuesta real a la pregunta está en su post (por favor, acredite su post en caso de duda). Estoy moviendo esta información anexa aquí para facilitar la recuperación de este post en el futuro. Sin restar importancia a ninguna de las otras excelentes respuestas, después de la respuesta de Mark, esto concluye la cuestión mediante la corrección del post en el OP.

Fuente

EN PITÓN:

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

2voto

El Six Puntos 89

Código Python:

import numpy as np

# desired correlation matrix
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)

# build some signals that will result in the desired correlation matrix
X = L.dot(np.random.normal(0,1,(3,1000))) # the more the sample (1000) the better

# estimate their correlation matrix
np.corrcoef(X)
array([[1.        , 0.58773667, 0.8978625 ],
       [0.58773667, 1.        , 0.47424997],
       [0.8978625 , 0.47424997, 1.        ]])

# Very good approxiamation :)

i-Ciencias.com

I-Ciencias es una comunidad de estudiantes y amantes de la ciencia en la que puedes resolver tus problemas y dudas.
Puedes consultar las preguntas de otros usuarios, hacer tus propias preguntas o resolver las de los demás.

Powered by:

X