Jak korzystać z rozkładu Cholesky'ego lub alternatywnie do skorelowanej symulacji danych

19

Używam dekompozycji Cholesky'ego do symulacji skorelowanych zmiennych losowych na podstawie macierzy korelacji. Chodzi o to, że wynik nigdy nie odtwarza struktury korelacji, jaka jest podana. Oto mały przykład w Pythonie ilustrujący sytuację.

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

To drukuje:

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

Jak widać, macierz korelacji szacowanej post hoc drastycznie różni się od poprzedniej. Czy w moim kodzie jest błąd, czy jest jakaś alternatywa dla użycia rozkładu Cholesky'ego?

Edytować

Proszę o wybaczenie za ten bałagan. Nie sądziłem, że wystąpił błąd w kodzie i / lub w sposobie, w jaki zastosowano rozkład Cholesky'ego z powodu jakiegoś nieporozumienia z materiałem, który studiowałem wcześniej. W rzeczywistości byłem pewien, że sama metoda nie była precyzyjna i nie przeszkadzało mi to aż do sytuacji, w której zadałem to pytanie. Dziękuję za wskazanie mojego błędnego przekonania. Zredagowałem tytuł, aby lepiej odzwierciedlić rzeczywistą sytuację zaproponowaną przez @Silverfish.

Eli Korvigo
źródło
1
Cholesky działa dobrze i jest to pytanie typu „czy możesz znaleźć błąd w moim kodzie”. Tytuł i treść pytania, jak zostało pierwotnie napisane, są w zasadzie „Cholesky nie działa, co jest alternatywą”? Będzie to bardzo mylące dla użytkowników przeszukujących tę stronę. Czy to pytanie należy edytować, aby to odzwierciedlić? (Minusem jest to, że odpowiedź javlacalle byłaby mniej istotna. Minusem jest to, że tekst pytania odzwierciedlałby wtedy to, co poszukiwacze faktycznie znaleźliby na stronie.)
Silverfish
@Antoni Parellada Tak, myślę, że przetłumaczyłeś mój kod MATLAB dla (a) poprawnego sposobu zrobienia tego w Python numpy, wraz z korektą dla np.linalg.cholesky będącego niższym trójkątem w porównaniu do cholery MATLABa jako trójkąta górnego. Przetłumaczyłem już niepoprawny kod OP na jego odpowiednik MATLAB i zduplikowałem jego nieprawidłowe wyniki.
Mark L. Stone,

Odpowiedzi:

11

Podejście oparte na rozkładzie Choleskiego powinno zadziałać, jest opisane tutaj i pokazane w odpowiedzi Mark L. Stone opublikowanej prawie w tym samym czasie, w którym ta odpowiedź.

N.(μ,Σ)

Y=QX+μ,zQ=Λ1/2)Φ,

YXΦΣΛΣΦ

Przykład w R(przepraszam, nie używam tego samego oprogramowania, którego użyłeś w pytaniu):

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

Możesz być także zainteresowany tym postem i tym postem .

javlacalle
źródło
Aby precyzyjnie odtworzyć macierz korelacji, przed zastosowaniem jej do procedury generowania danych należy usunąć fałszywe korelacje w danych losowych z generatora losowego. Na przykład sprawdź korelację danych losowych w eps, aby najpierw zobaczyć te fałszywe korelacje.
Gottfried Helms
17

Ludzie mogliby znaleźć błąd znacznie szybciej, gdybyś wyjaśnił, co zrobiłeś słowami i algebrą, a nie kodem (lub przynajmniej pisząc go za pomocą pseudokodu).

Wygląda na to, że robisz odpowiednik tego (choć prawdopodobnie transponowanego):

  1. n×kZ

  2. σjaμja

  3. Y=L.X

L.

Co powinieneś zrobić, to:

  1. n×kZ

  2. X=L.Z

  3. σjaμja

Istnieje wiele wyjaśnień tego algorytmu na stronie. na przykład

Jak wygenerować skorelowane liczby losowe (podane średnie, wariancje i stopień korelacji)?

Czy mogę użyć metody Cholesky'ego do generowania skorelowanych zmiennych losowych o podanej średniej?

Ten omawia go bezpośrednio pod kątem pożądanej macierzy kowariancji, a także podaje algorytm uzyskiwania pożądanej kowariancji próbki :

Generowanie danych za pomocą danej macierzy kowariancji próbki

Glen_b - Przywróć Monikę
źródło
11

Nie ma nic złego w faktoryzacji Cholesky'ego. W twoim kodzie jest błąd. Zobacz edycję poniżej.

Oto kod MATLAB i wyniki, najpierw dla n_obs = 10000, jak masz, a następnie dla n_obs = 1e8. Dla uproszczenia, ponieważ nie wpływa to na wyniki, nie zawracam sobie głowy środkami, tzn. Robię je zerami. Zauważ, że chol MATLABa wytwarza górny trójkątny czynnik Cholesky'ego R matrycy M, tak że R '* R = M. numpy.linalg.cholesky wytwarza niższy trójkątny czynnik Cholesky'ego, więc potrzebna jest korekta względem mojego kodu; ale uważam, że pod tym względem twój kod jest w porządku.

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

Edycja: Znalazłem twój błąd. Niepoprawnie zastosowałeś odchylenie standardowe. Jest to odpowiednik tego, co zrobiłeś, co jest złe.

   >> 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
Mark L. Stone
źródło
6

W CV nie chodzi o kod, ale byłem zaintrygowany, aby zobaczyć, jak to będzie wyglądać po wszystkich dobrych odpowiedziach, a konkretnie w @Mark L. Stone. Rzeczywista odpowiedź na to pytanie znajduje się na jego stanowisku (w razie wątpliwości prosimy o uznanie go) Przenoszę tutaj te dołączone informacje, aby ułatwić pobieranie tego postu w przyszłości. Bez pomijania innych doskonałych odpowiedzi, po odpowiedzi Marka, kończy to problem, poprawiając post w PO.

Źródło

W 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
Antoni Parellada
źródło
1

Jak inni już pokazali: cholesky działa. Oto fragment kodu, który jest bardzo krótki i bardzo zbliżony do pseudokodu: codepiece w 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
Gottfried Helms
źródło