Jak efektywnie generować losowe macierze korelacji dodatnio-półpokrytych?

38

Chciałbym być w stanie efektywnie generować macierze korelacji dodatnich-półprzewodnikowych (PSD). Moja metoda gwałtownie zwalnia, gdy zwiększam rozmiar generowanych macierzy.

  1. Czy możesz zasugerować jakieś skuteczne rozwiązania? Jeśli znasz jakieś przykłady w Matlabie, byłbym bardzo wdzięczny.
  2. W jaki sposób przy generowaniu macierzy korelacji PSD wybrałbyś parametry opisujące macierze do wygenerowania? Średnia korelacja, standardowe odchylenie korelacji, wartości własne?
Eduardas
źródło

Odpowiedzi:

16

Możesz to zrobić wstecz: każda macierz (zbiór wszystkich symetrycznych macierzy PSD) może zostać rozłożona jako p × pCR++pp×p

C=OTDO gdzie jest macierzą ortonormalnąO

Aby uzyskać , najpierw wygeneruj losową podstawę (gdzie są wektorami losowymi, zwykle w ). Stamtąd użyj procesu ortogonalizacji Gram-Schmidta, aby uzyskać( V 1 , . . . , V s ) v I ( - 1 , 1 ) ( U 1 , . . . . , U t ) = OO(v1,...,vp)vi(1,1)(u1,....,up)=O

R ma wiele pakietów, które mogą efektywnie wykonywać losową ortogonalizację GS, to znaczy nawet dla dużych wymiarów, na przykład pakiet „daleki”. Chociaż algorytm GS znajdziesz na wiki, prawdopodobnie lepiej nie wymyślać na nowo koła i wybierać implementację Matlaba (jedna z pewnością istnieje, po prostu nie mogę jej polecić).

Wreszcie, jest macierzami ukośnymi, których elementy są dodatnie (znowu jest to łatwe do wygenerowania: wygeneruj liczb losowych, wyprostuj je, posortuj i umieść na przekątnej macierzy po ).p p pDppp

użytkownik603
źródło
3
(1) Należy zauważyć, że wynikowe nie będzie macierzą korelacji (zgodnie z żądaniem PO), ponieważ nie będzie miało macierzy na przekątnej. Oczywiście mogą być przeskalowane mieć te po przekątnej, przez ustawienie go na , gdzie jest macierzą diagonalną o tej samej przekątnej, jak . (2) Jeśli się nie mylę, spowoduje to powstanie macierzy korelacji, w których wszystkie elementy o przekątnej są skoncentrowane wokół , więc nie ma elastyczności, której szukał PO (OP chciał mieć możliwość ustawienia „średniej korelacji” , standardowe odchylenie korelacji, wartości własne ” )E - 1 / 2 C, E - 1 / 2 e C 0CE1/2CE1/2EC0
mówi ameba Przywróć Monikę
@amoeba: Zajmę się (2), ponieważ, jak zauważyłeś, rozwiązanie (1) jest trywialne. Jedną liczbą charakteryzującą „kształt” (związek między wejściem i wyjściem elementów diagonalnych) macierzy PSD (a zatem kowariancją, a także macierzą korelacji) jest jej liczba warunkowa. I metoda opisana powyżej pozwala na pełną kontrolę nad tym. „Stężenie elementów niesąsych wokół 0” nie jest cechą metody stosowanej do generowania matryc PSD, ale raczej konsekwencją ograniczeń koniecznych do zapewnienia, że ​​macierz jest PSD i fakt, że jest duże. p
user603,
Czy mówisz, że wszystkie duże matryce PSD mają elementy o przekątnej zbliżonej do zera? Nie zgadzam się, tak nie jest. Sprawdź moją odpowiedź tutaj, aby znaleźć kilka przykładów: Jak wygenerować macierz losowej korelacji, która ma w przybliżeniu normalnie rozmieszczone wpisy poziagonalne przy danym odchyleniu standardowym? Ale można bezpośrednio zobaczyć, że nie jest to przypadek, ponieważ macierzą kwadratową posiadające wszystkie jedynki na przekątnej i stałej wartości wszędzie poza przekątna jest PSD i może być dowolnie duża (ale oczywiście poniżej ). ρ 1ρρ1
ameba mówi Przywróć Monikę
@amoeba: wtedy pomyliłem się, zakładając, że z konieczności nierównomierna duża macierz korelacji, gdy mogą być zarówno dodatnie, jak i ujemne, jest bliska 0. Dzięki za oświecający przykład.
user603,
1
Przeczytałem bardzo fajny artykuł na temat generowania macierzy losowej korelacji i podałem tutaj swoją własną odpowiedź (a także inną odpowiedź w tym powiązanym wątku). Myślę, że może Cię to zainteresować.
ameba mówi Przywróć Monikę
27

Artykuł Generowanie losowych macierzy korelacji na podstawie metody winorośli i rozszerzonej cebuli autorstwa Lewandowskiego, Kurowickiej i Joe (LKJ), 2009, zapewnia ujednolicone traktowanie i opis dwóch skutecznych metod generowania macierzy losowej korelacji. Obie metody pozwalają na generowanie macierzy z jednolitego rozkładu w ściśle określonym sensie zdefiniowanym poniżej, są łatwe do wdrożenia, szybkie i mają dodatkową zaletę polegającą na zabawnych nazwach.

Prawdziwa macierz symetryczna o rozmiarze z elementami na przekątnej ma unikalne elementy poza przekątną, więc można ją sparametryzować jako punkt w . Każdy punkt w tej przestrzeni odpowiada macierzy symetrycznej, ale nie wszystkie z nich są dodatnio określone (ponieważ muszą to być macierze korelacji). Macierze korelacji tworzą zatem podzbiór (właściwie połączony wypukły podzbiór), i obie metody mogą generować punkty z jednolitego rozkładu w tym podzbiorze.d×dd(d1)/2Rd(d1)/2Rd(d1)/2

Przedstawię własną implementację MATLAB dla każdej metody i zilustruję je .d=100


Metoda cebuli

Metoda cebulowa pochodzi z innego artykułu (nr 3 w LKJ) i zawdzięcza swoją nazwę faktowi, że macierze korelacji są generowane zaczynając od macierzy i rosnąc ją kolumna po kolumnie i rząd po rzędzie. Wynikowy rozkład jest jednolity. Naprawdę nie rozumiem matematyki za tą metodą (i tak wolę drugą metodę), ale oto wynik:1×1

Metoda cebuli

Tu i poniżej tytuł każdej podploty pokazuje najmniejsze i największe wartości własne oraz wyznacznik (iloczyn wszystkich wartości własnych). Oto kod:

%// ONION METHOD to generate random correlation matrices distributed randomly
function S = onion(d)
    S = 1;
    for k = 2:d
        y = betarnd((k-1)/2, (d-k)/2); %// sampling from beta distribution
        r = sqrt(y);
        theta = randn(k-1,1);
        theta = theta/norm(theta);
        w = r*theta;
        [U,E] = eig(S);
        R = U*E.^(1/2)*U';             %// R is a square root of S
        q = R*w;
        S = [S q; q' 1];               %// increasing the matrix size
    end
end

Metoda rozszerzonej cebuli

LKJ nieznacznie zmodyfikuje tę metodę, aby móc próbkować macierze korelacji z rozkładu proporcjonalnego do . Im większa , tym większa będzie wyznacznik, co oznacza, że ​​generowane macierze korelacji będą coraz bardziej zbliżały się do matrycy tożsamości. Wartość odpowiada rozkładowi jednorodnemu. Na poniższym rysunku macierze są generowane z .C[detC]η1ηη=1η=1,10,100,1000,10000,100000

Metoda rozszerzonej cebuli

Z jakiegoś powodu, aby uzyskać wyznacznik tego samego rzędu wielkości co w metodzie z wanilią, muszę a nie (jak twierdzi LKJ). Nie jestem pewien, gdzie jest błąd.η=0η=1

%// EXTENDED ONION METHOD to generate random correlation matrices
%// distributed ~ det(S)^eta [or maybe det(S)^(eta-1), not sure]
function S = extendedOnion(d, eta)
    beta = eta + (d-2)/2;
    u = betarnd(beta, beta);
    r12 = 2*u - 1;
    S = [1 r12; r12 1];  

    for k = 3:d
        beta = beta - 1/2;
        y = betarnd((k-1)/2, beta);
        r = sqrt(y);
        theta = randn(k-1,1);
        theta = theta/norm(theta);
        w = r*theta;
        [U,E] = eig(S);
        R = U*E.^(1/2)*U';
        q = R*w;
        S = [S q; q' 1];
    end
end

Metoda winorośli

Metoda Vine została pierwotnie zasugerowana przez Joe (J w LKJ) i ulepszona przez LKJ. Bardziej mi się podoba, ponieważ jest koncepcyjnie łatwiejszy, a także łatwiejszy do modyfikacji. Chodzi o to, aby wygenerować częściowe korelacje (są one niezależne i mogą mieć dowolne wartości z bez żadnych ograniczeń), a następnie przekształcić je w surowe korelacje za pomocą wzoru rekurencyjnego. Wygodnie jest zorganizować obliczenia w określonej kolejności, a ten wykres jest znany jako „winorośl”. Co ważne, jeśli próbkuje się korelacje częściowe z poszczególnych rozkładów beta (różnych dla różnych komórek w macierzy), wówczas uzyskana macierz będzie rozkładana równomiernie. Tutaj ponownie LKJ wprowadza dodatkowy parametr aby pobrać próbkę z rozkładu proporcjonalnego dod(d1)/2[1,1]η[detC]η1 . Wynik jest identyczny z rozszerzoną cebulą:

Metoda winorośli

%// VINE METHOD to generate random correlation matrices
%// distributed ~ det(S)^eta [or maybe det(S)^(eta-1), not sure]
function S = vine(d, eta)
    beta = eta + (d-1)/2;   
    P = zeros(d);           %// storing partial correlations
    S = eye(d);

    for k = 1:d-1
        beta = beta - 1/2;
        for i = k+1:d
            P(k,i) = betarnd(beta,beta); %// sampling from beta
            P(k,i) = (P(k,i)-0.5)*2;     %// linearly shifting to [-1, 1]
            p = P(k,i);
            for l = (k-1):-1:1 %// converting partial correlation to raw correlation
                p = p * sqrt((1-P(l,i)^2)*(1-P(l,k)^2)) + P(l,i)*P(l,k);
            end
            S(k,i) = p;
            S(i,k) = p;
        end
    end
end

Metoda winorośli z ręcznym próbkowaniem korelacji cząstkowych

Jak widać powyżej, równomierny rozkład daje prawie diagonalne macierze korelacji. Ale można łatwo zmodyfikować metodę winorośli, aby uzyskać silniejsze korelacje (nie jest to opisane w pracy LKJ, ale jest proste): w tym celu należy próbkować korelacje częściowe z rozkładu skoncentrowanego wokół . Poniżej próbuję je z rozkładu beta (przeskalowanego z na ) za pomocą . Im mniejsze parametry rozkładu beta, tym bardziej jest ono skoncentrowane w pobliżu krawędzi.±1[0,1][1,1]α=β=50,20,10,5,2,1

Metoda winorośli z ręcznym pobieraniem próbek

Zauważ, że w tym przypadku nie ma gwarancji, że rozkład będzie niezmienny w permutacji, więc dodatkowo generuję losowo wiersze i kolumny po wygenerowaniu.

%// VINE METHOD to generate random correlation matrices
%// with all partial correlations distributed ~ beta(betaparam,betaparam)
%// rescaled to [-1, 1]
function S = vineBeta(d, betaparam)
    P = zeros(d);           %// storing partial correlations
    S = eye(d);

    for k = 1:d-1
        for i = k+1:d
            P(k,i) = betarnd(betaparam,betaparam); %// sampling from beta
            P(k,i) = (P(k,i)-0.5)*2;     %// linearly shifting to [-1, 1]
            p = P(k,i);
            for l = (k-1):-1:1 %// converting partial correlation to raw correlation
                p = p * sqrt((1-P(l,i)^2)*(1-P(l,k)^2)) + P(l,i)*P(l,k);
            end
            S(k,i) = p;
            S(i,k) = p;
        end
    end

    %// permuting the variables to make the distribution permutation-invariant
    permutation = randperm(d);
    S = S(permutation, permutation);
end

Oto, jak histogramy elementów nie przekątnych wyglądają dla powyższych matryc (wariancja rozkładu monotonicznie wzrasta):

Elementy o przekątnej


Aktualizacja: przy użyciu losowych czynników

Jedna bardzo prosta metoda generowania macierzy losowych korelacji z pewnymi silnymi korelacjami została zastosowana w odpowiedzi przez @shabbychef i chciałbym to również zilustrować tutaj. Chodzi o to, aby losowo wygenerować kilka ( ) ładunków czynnikowych (losowa macierz wielkości), utworzyć macierz kowariancji (który oczywiście nie będzie pełnej rangi ) i dodaj do niej losową macierz ukośną z dodatnimi elementami, aby pełna ranga. Otrzymaną macierz kowariancji można znormalizować, aby stała się macierzą korelacji, pozwalającW k x d W WD B = W W + D C = E - 1 / 2 B E - 1 / 2 PL B K = 100 , 50 , 20 , 10 , 5 , 1k<dWk×dWWDB=WW+DC=E1/2BE1/2Gdzie jest macierzą diagonalną o tej samej przekątnej, co . To jest bardzo proste i załatwia sprawę. Oto kilka przykładowych macierzy korelacji dla :EBk=100,50,20,10,5,1

losowe macierze korelacji z czynników losowych

I kod:

%// FACTOR method
function S = factor(d,k)
    W = randn(d,k);
    S = W*W' + diag(rand(1,d));
    S = diag(1./sqrt(diag(S))) * S * diag(1./sqrt(diag(S)));
end

Oto kod opakowania używany do generowania liczb:

d = 100; %// size of the correlation matrix

figure('Position', [100 100 1100 600])
for repetition = 1:6
    S = onion(d);

    %// etas = [1 10 100 1000 1e+4 1e+5];
    %// S = extendedOnion(d, etas(repetition));

    %// S = vine(d, etas(repetition));

    %// betaparams = [50 20 10 5 2 1];
    %// S = vineBeta(d, betaparams(repetition));

    subplot(2,3,repetition)

    %// use this to plot colormaps of S
    imagesc(S, [-1 1])
    axis square
    title(['Eigs: ' num2str(min(eig(S)),2) '...' num2str(max(eig(S)),2) ', det=' num2str(det(S),2)])

    %// use this to plot histograms of the off-diagonal elements
    %// offd = S(logical(ones(size(S))-eye(size(S))));
    %// hist(offd)
    %// xlim([-1 1])
end
ameba mówi Przywróć Monikę
źródło
2
To fantastyczne podsumowanie. Cieszę się, że coś powiedziałem!
shadowtalker
Kiedy przetłumaczyłem kod matlab macierzy korelacji opartej na winorośli na R i przetestowałem go, gęstość korelacji w kolumnie 1 była zawsze inna niż w późniejszych kolumnach. Być może przetłumaczyłem coś niepoprawnie, ale być może ta notatka pomaga komuś.
Charlie,
3
Dla użytkowników R funkcja rcorrmatrix w pakiecie klasterGeneration (napisana przez W Qui i H. Joe) implementuje metodę winorośli.
RNM,
15

Jeszcze prostszym charakterystyki jest to, że dla rzeczywistej macierzy , jest dodatnia półokreśloną. Aby zrozumieć, dlaczego tak jest, wystarczy udowodnić, że dla wszystkich wektorów (oczywiście o odpowiednim rozmiarze). Jest to trywialne:co jest nieujemne. Więc w Matlabie po prostu spróbujAATAyT(ATA)y0yyT(ATA)y=(Ay)TAy=||Ay||

A = randn(m,n);   %here n is the desired size of the final matrix, and m > n
X = A' * A;

W zależności od aplikacji może to nie zapewniać odpowiedniej dystrybucji wartości własnych; Odpowiedź Kwaka jest pod tym względem znacznie lepsza. Wartości własne Xwytworzone przez ten fragment kodu powinny być zgodne z rozkładem Marchenko-Pastur.

Powiedzmy, że do symulacji macierzy korelacji zapasów może być potrzebne nieco inne podejście:

k = 7;      % # of latent dimensions;
n = 100;    % # of stocks;
A = 0.01 * randn(k,n);  % 'hedgeable risk'
D = diag(0.001 * randn(n,1));   % 'idiosyncratic risk'
X = A'*A + D;
ascii_hist(eig(X));    % this is my own function, you do a hist(eig(X));
-Inf <= x <  -0.001 : **************** (17)
-0.001 <= x <   0.001 : ************************************************** (53)
 0.001 <= x <   0.002 : ******************** (21)
 0.002 <= x <   0.004 : ** (2)
 0.004 <= x <   0.005 :  (0)
 0.005 <= x <   0.007 : * (1)
 0.007 <= x <   0.008 : * (1)
 0.008 <= x <   0.009 : *** (3)
 0.009 <= x <   0.011 : * (1)
 0.011 <= x <     Inf : * (1)
shabbychef
źródło
1
czy byłbyś skłonny podzielić się swoją funkcją ascii_hist przez przypadek?
btown
@btown margines jest zbyt mały, aby go pomieścić!
shabbychef
1
Literówka występuje w- brakuje jej ostatniego kwadratu! yT(ATA)y=(Ay)TAy=||Ay||
Silverfish,
8

W odmianie odpowiedzi kwaka: wygeneruj macierz diagonalną z losowymi nieujemnymi wartościami własnymi z wybranego rozkładu, a następnie wykonaj transformację podobieństwa z pseudolosową macierzą rozproszoną przez Haara .A = Q D Q T QDA=QDQTQ

JM nie jest statystykiem
źródło
M.: Niezłe odniesienie: wydaje się, że jest to najbardziej wydajne rozwiązanie (asymptotycznie).
whuber
3
@whuber: Heh, odebrałem go od Goluba i Van Loana (oczywiście); Używam tego cały czas, aby pomóc w generowaniu macierzy testowych do testów warunków własnych wartości własnych / procedur wartości pojedynczych. Jak widać z pracy, jest to w zasadzie odpowiednik rozkładania QR przypadkowej macierzy, jak sugerował kwak, z tym wyjątkiem, że odbywa się to bardziej wydajnie. Implementacja MATLAB jest dostępna w Higham's Text Matrix Toolbox, BTW.
JM nie jest statystykiem
M .:> Dziękujemy za wdrożenie Matlaba. Czy przypadkiem znasz pseudolosowy generator macierzy Haar w języku R?
user603
@kwak: Nie mam pojęcia, ale jeśli nie ma jeszcze implementacji, przetłumaczenie kodu MATLAB na R nie powinno być zbyt trudne (mogę spróbować go poderwać, jeśli naprawdę go nie ma); jedynym warunkiem jest przyzwoity generator normalnych zmiennych pseudolosowych, który jestem pewien, że R. ma.
JM nie jest statystykiem
M .:> tak, prawdopodobnie przetłumaczę to sobie. Dzięki za linki, Best.
user603
4

Nie określono rozkładu macierzy. Dwa popularne to Wishart i odwrotne rozkłady Wishart. Rozkładu Bartlett daje Cholesky'iego faktoryzacji macierzy losowo Wishart (która może być również skutecznie rozwiązać w celu uzyskania losowo odwrotną Wishart matrycy).

W rzeczywistości przestrzeń Choleskiego jest wygodnym sposobem generowania innych rodzajów losowych macierzy PSD, ponieważ wystarczy upewnić się, że przekątna nie jest ujemna.

Simon Byrne
źródło
> Nieprzypadkowo: dwie macierze wygenerowane z tego samego Whisharda nie będą od siebie niezależne. Jeśli planujesz zmienić Whishart w każdym pokoleniu, to jak zamierzasz generować Whishart w pierwszej kolejności?
user603 24.09.10
@kwak: Nie rozumiem twojego pytania: rozkład Bartletta da niezależne losowania z tej samej dystrybucji Wishart.
Simon Byrne
> Pozwólcie, że sformułuję to, skąd czerpiesz macierz skali swojej dystrybucji whishart?
user603,
1
@kwak: jest to parametr rozkładu, więc jest ustalony. Wybierasz go na początku, na podstawie pożądanych cech twojego rozkładu (takich jak średnia).
Simon Byrne
3

Najprostszą metodą jest powyższa, która jest symulacją losowego zestawu danych i obliczeniem Gramiana . Słowo ostrzeżenia: wynikowa macierz nie będzie jednorodnie losowa, ponieważ jej rozkład, powiedzmy, będzie miał obroty nie rozmieszczone zgodnie z miarą Haara. Jeśli chcesz mieć „równomiernie rozmieszczone” matryce PSD, możesz zastosować dowolne z opisanych tutaj podejść .UTSU

niezadowolony
źródło
Jeśli wpisy są generowane z rozkładu normalnego, a nie z jednorodnego, wspomniany rozkład powinien być niezmiennikiem SO (n) (a zatem równomiernie rozłożony względem miary Haara).
whuber
Ciekawy. Czy możesz podać referencje?
gappy
1
> Problem z tą metodą polega na tym, że nie można kontrolować stosunku najmniejszej do największej wartości własnej (i myślę, że ponieważ rozmiar losowo wygenerowanego zestawu danych osiągnie nieskończoność, stosunek ten zbiegnie się do 1).
user603,
1

Jeśli chcesz mieć większą kontrolę nad wygenerowaną symetryczną matrycą PSD, np. Wygenerować zestaw danych do weryfikacji syntetycznej, masz do dyspozycji szereg parametrów. Symetryczna macierz PSD odpowiada hiper-elipsie w przestrzeni N-wymiarowej, ze wszystkimi powiązanymi stopniami swobody:

  1. Rotacje
  2. Długości osi.

Tak więc dla macierzy dwuwymiarowej (tj. Elipsy 2d) będziesz mieć 1 obrót + 2 osie = 3 parametry.

Σ=ODOTΣOD

Σ

figure;
mu = [0,0];
for i=1:16
    subplot(4,4,i)
    theta = (i/16)*2*pi;   % theta = rand*2*pi;
    U=[cos(theta), -sin(theta); sin(theta) cos(theta)];
    % The diagonal's elements control the lengths of the axes
    D = [10, 0; 0, 1]; % D = diag(rand(2,1));    
    sigma = U*D*U';
    data = mvnrnd(mu,sigma,1000);
    plot(data(:,1),data(:,2),'+'); axis([-6 6 -6 6]); hold on;
end

U

PeriRamm
źródło
0

Tanim i wesołym podejściem, którego użyłem do testowania, jest wygenerowanie m N (0,1) n wektorów V [k], a następnie użycie P = d * I + Suma {V [k] * V [k] '} jako macierz psd nxn. Przy m <n będzie to liczba pojedyncza dla d = 0, a dla małego d będzie miała wysoką liczbę warunkową.


źródło
2
> Problem z tą metodą polega na tym, że nie można kontrolować stosunku najmniejszej do największej wartości własnej (i myślę, że ponieważ rozmiar losowo wygenerowanego zestawu danych osiągnie nieskończoność, stosunek ten zbiegnie się do 1).
user603,
> poza tym metoda nie jest bardzo wydajna (z obliczeniowego punktu widzenia)
user603 24.09.10
1
Twoja „macierz losowa” jest specjalnie skonstruowana i nazywana „macierzą diagonalną plus ranga 1” (macierz DR1), więc nie jest tak naprawdę dobrą reprezentatywną macierzą losową.
JM nie jest statystykiem