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

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?

Możesz to zrobić wstecz: każda macierz (zbiór wszystkich symetrycznych macierzy PSD) może zostać rozłożona jako p × $C \in \mathbb{R}_{++}^p$$p \times p$

$C=O^{T}DO$ 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 ) = $O$$(v_1,...,v_p)$$v_i$$(-1,1)$$(u_1,....,u_p)=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 $D$$p$$p$$p$

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 \times d$$d(d-1)/2$$\mathbb R^{d(d-1)/2}$$\mathbb R^{d(d-1)/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\times 1$

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 .$\mathbf C$$[\mathrm{det}\:\mathbf C]^{\eta-1}$$\eta$$\eta=1$$\eta={1, 10, 100, 1000, 10\:000, 100\:000}$

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.$\eta=0$$\eta=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 do$d(d-1)/2$$[-1, 1]$$\eta$$[\mathrm{det}\:\mathbf C]^{\eta-1}$ . Wynik jest identyczny z rozszerzoną cebulą:

%// 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.$\pm 1$$[0,1]$$[-1, 1]$$\alpha=\beta={50, 20, 10, 5, 2, 1}$

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

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 W ⊤ D B = W W ⊤ + D C = E - 1 / 2 B E - 1 / 2 PL B K = 100 , 50 , 20 , 10 , 5 , $k<d$$\mathbf W$$k \times d$$\mathbf W \mathbf W^\top$$\mathbf D$$\mathbf B = \mathbf W \mathbf W^\top + \mathbf D$$\mathbf C = \mathbf E^{-1/2}\mathbf B \mathbf E^{-1/2}$Gdzie jest macierzą diagonalną o tej samej przekątnej, co . To jest bardzo proste i załatwia sprawę. Oto kilka przykładowych macierzy korelacji dla :$\mathbf E$$\mathbf B$$k={100, 50, 20, 10, 5, 1}$

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

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óbuj$A$$A^TA$$y^T (A^TA) y \ge 0$$y$$y^T (A^TA) y = (Ay)^T Ay = ||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)

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 $\mathbf{D}$$\mathbf{A}=\mathbf{Q}\mathbf{D}\mathbf{Q}^T$$\mathbf{Q}$

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.

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ść .$U^TSU$

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.

$\Sigma=ODO^T$$\Sigma$$O$$D$

$\Sigma$

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$

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