Verwendung der Cholesky-Zerlegung oder einer Alternative für die Simulation korrelierter Daten

Ich benutze die Cholesky-Zerlegung, um korrelierte Zufallsvariablen bei gegebener Korrelationsmatrix zu simulieren. Die Sache ist, das Ergebnis reproduziert niemals die Korrelationsstruktur, wie sie gegeben ist. Hier ist ein kleines Beispiel in Python, um die Situation zu veranschaulichen.

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

Dies druckt:

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

Wie Sie sehen, unterscheidet sich die post-hoc geschätzte Korrelationsmatrix drastisch von der vorherigen. Befindet sich in meinem Code ein Fehler oder gibt es eine Alternative zur Verwendung der Cholesky-Dekomposition?

Bearbeiten

Ich bitte um Verzeihung für dieses Durcheinander. Ich dachte nicht, dass es einen Fehler im Code und / oder in der Art und Weise gab, wie die Cholesky-Zerlegung angewendet wurde, weil ich das Material, das ich zuvor studiert hatte, falsch verstanden hatte. Tatsächlich war ich mir sicher, dass die Methode selbst nicht genau sein sollte, und ich war damit bis zu der Situation einverstanden, die mich veranlasste, diese Frage zu stellen. Vielen Dank für den Hinweis auf das Missverständnis, das ich hatte. Ich habe den Titel bearbeitet, um die von @Silverfish vorgeschlagene reale Situation besser widerzuspiegeln.

Der Ansatz, der auf der Cholesky-Zerlegung basiert, sollte funktionieren, er wird hier beschrieben und in der Antwort von Mark L. Stone fast zeitgleich mit dieser Antwort veröffentlicht.

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

wo die letzten Züge sind, $Y$$X$$\Phi$$\Sigma$$\Lambda$$\Sigma$$\Phi$

Beispiel in R(Entschuldigung, ich verwende nicht dieselbe Software, die Sie in der Frage verwendet haben):

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

Sie könnten auch an diesem Beitrag und diesem Beitrag interessiert sein .

Die Leute würden Ihren Fehler wahrscheinlich viel schneller finden, wenn Sie erklären würden, was Sie mit Worten und Algebra anstatt mit Code gemacht haben (oder zumindest mit Pseudocode schreiben).

Sie scheinen das Äquivalent dazu zu tun (obwohl möglicherweise transponiert):

1. Generiere ein $n\times k$$Z$

2. multiplizieren Sie die Spalten mit $\sigma_i$$\mu_i$

3. berechne $Y=LX$

$L$

Was Sie tun sollten, ist Folgendes:

1. Generiere ein $n\times k$$Z$

2. berechne $X=LZ$

3. multipliziere die Spalten mit und addiere$\sigma_i$$\mu_i$

Es gibt viele Erklärungen für diesen Algorithmus vor Ort. z.B

Wie werden korrelierte Zufallszahlen generiert (gegebene Mittelwerte, Varianzen und Grad der Korrelation)?

Kann ich die Cholesky-Methode verwenden, um korrelierte Zufallsvariablen mit dem angegebenen Mittelwert zu erzeugen?

Dieser diskutiert es direkt in Bezug auf die gewünschte Kovarianzmatrix und gibt auch einen Algorithmus zum Erhalten einer gewünschten Sample- Kovarianz an:

Generieren von Daten mit einer bestimmten Stichproben-Kovarianzmatrix

An der Cholesky-Faktorisierung ist nichts auszusetzen. In Ihrem Code ist ein Fehler aufgetreten. Siehe unten bearbeiten.

Hier ist der MATLAB-Code und die Ergebnisse, zuerst für n_obs = 10000 wie Sie, dann für n_obs = 1e8. Der Einfachheit halber kümmere ich mich nicht um die Mittelwerte, da sie die Ergebnisse nicht beeinflussen, dh ich mache sie zu Nullen. Beachten Sie, dass MATLABs chol einen oberen dreieckigen Cholesky-Faktor R der Matrix M erzeugt, so dass R '* R = M. numpy.linalg.cholesky einen unteren dreieckigen Cholesky-Faktor erzeugt, so dass eine Anpassung gegenüber meinem Code erforderlich ist; aber ich glaube, Ihr Code ist in dieser Hinsicht in Ordnung.

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

Edit: Ich habe deinen Fehler gefunden. Sie haben die Standardabweichung falsch angewendet. Dies ist das Äquivalent zu dem, was Sie getan haben, was falsch ist.

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

In meinem Lebenslauf geht es nicht um Code, aber ich war gespannt, wie dies nach all den guten Antworten und speziell nach dem Beitrag von @Mark L. Stone aussehen würde. Die eigentliche Antwort auf die Frage finden Sie auf seinem Beitrag (im Zweifelsfall bitte seinen Beitrag gutschreiben). Ich verschiebe diese angehängten Informationen hierher, um den Abruf dieses Beitrags in Zukunft zu erleichtern. Ohne eine der anderen hervorragenden Antworten herunterzuspielen, wird das Problem nach Marks Antwort durch die Korrektur des Postens im OP behoben.

Quelle

IN 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.        ]]

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

Wie schon andere gezeigt haben: cholesky works. Hier ein Teil des Codes, der sehr kurz und sehr nahe am Pseudocode ist: ein Codeteil in 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