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.
źródło
Odpowiedzi:
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ź.
Przykład w
R
(przepraszam, nie używam tego samego oprogramowania, którego użyłeś w pytaniu):Możesz być także zainteresowany tym postem i tym postem .
źródło
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):
Co powinieneś zrobić, to:
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
źródło
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.
Edycja: Znalazłem twój błąd. Niepoprawnie zastosowałeś odchylenie standardowe. Jest to odpowiednik tego, co zrobiłeś, co jest złe.
źródło
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:
IN [R]:
źródło
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:
źródło