Jak wygenerować dużą macierz korelacji losowej pełnej rangi z obecnymi silnymi korelacjami?

Chciałbym wygenerować losową macierz korelacji o wielkości n × n, tak aby występowały pewne umiarkowanie silne korelacje:$\mathbf C$$n \times n$

• kwadratowa rzeczywista macierz symetryczna o rozmiarze , np. n = 100 ;$n \times n$$n=100$
• pozytywnie określone, tj. ze wszystkimi wartościami własnymi rzeczywistymi i dodatnimi;
• pełna ranga;
• wszystkie elementy ukośne równe ;$1$
• elementy o przekątnej powinny być odpowiednio równomiernie rozmieszczone na . Dokładny rozkład nie ma znaczenia, ale chciałbym mieć jakąś umiarkowanie dużą ilość (np. 10 % ) umiarkowanie dużych wartości (np. O wartości bezwzględnej 0,5 lub wyższej). Zasadniczo chcę się upewnić, że C nie jest prawie przekątna ze wszystkimi elementami nie przekątnymi ≈ 0 .$(-1, 1)$$10\%$$0.5$$\mathbf C$$\approx 0$

Czy istnieje prosty sposób to zrobić?

Celem jest wykorzystanie takich losowych macierzy do porównania niektórych algorytmów pracujących z macierzami korelacji (lub kowariancji).

Metody, które nie działają

Oto kilka sposobów generowania losowych macierzy korelacji, o których wiem, ale które tutaj nie działają:

1. Wygeneruj losowy o rozmiarze s × n , wyśrodkuj, standaryzuj i utwórz macierz korelacji C = $\mathbf X$$s \times n$. Jeślis>n, spowoduje to na ogół, że wszystkie korelacje poza przekątną będą wynosić około0. Jeślis≪n, niektóre korelacje będą silne, aleCnie będzie pełnej rangi.$\mathbf C=\frac{1}{s-1}\mathbf X^\top \mathbf X$$s>n$$0$$s\ll n$$\mathbf C$

2. Wygeneruj losową dodatnią określoną macierz na jeden z następujących sposobów:$\mathbf B$

   • Wygeneruj losowy kwadrat i zrób symetryczny dodatni określony B = A A ⊤ .$\mathbf A$$\mathbf B = \mathbf A \mathbf A^\top$

   • Wygeneruj losowy kwadrat , zrób symetryczny E = A + A ⊤ i ustaw go jako dodatni, wykonując rozkład własny E = U S U ⊤ i ustawiając wszystkie ujemne wartości własne na zero: B = $\mathbf A$$\mathbf E = \mathbf A + \mathbf A^\top$$\mathbf E = \mathbf U \mathbf S \mathbf U^\top$ . Uwaga: spowoduje to macierz z niedoborem rang.$\mathbf B = \mathbf U \:\mathrm{max}\{\mathbf S, \mathbf 0\} \:\mathbf U^\top$

   • Generuj losowy ortogonalny (np. Generując losowy kwadrat A i dokonując jego rozkładu QR lub w procesie Grama-Schmidta) i losową przekątną D ze wszystkimi dodatnimi elementami; forma B = Q D Q ⊤ .$\mathbf Q$$\mathbf A$$\mathbf D$$\mathbf B = \mathbf Q \mathbf D \mathbf Q^\top$

   Otrzymany macierz może być łatwo znormalizowane, aby wszystkie te diagonalne: C = D - 1 / 2 B D - 1 / 2 , gdzie D = d i $\mathbf B$$\mathbf C = \mathbf D^{-1/2}\mathbf B \mathbf D^{-1/2}$ jest macierzą diagonalną o tej samej przekątnej jako B . Wszystkie trzy wymienione powyżej sposoby generowania B powodują,że C ma elementy o przekątnej zamkniętej 0 .$\mathbf D = \mathrm{diag}\:\mathbf B$$\mathbf B$$\mathbf B$$\mathbf C$$0$

Aktualizacja: starsze wątki

Po opublikowaniu pytania znalazłem w przeszłości dwa prawie duplikaty:

• Jak wygenerować macierz losowej korelacji, która ma w przybliżeniu normalnie rozmieszczone wpisy o przekątnej przy danym odchyleniu standardowym?
• Jak efektywnie generować losowe macierze korelacji dodatnio-półpokrytych?

Niestety żaden z tych wątków nie zawiera satysfakcjonującej odpowiedzi (do tej pory :)

Inne odpowiedzi wymyśliły fajne sztuczki, aby rozwiązać mój problem na różne sposoby. Jednak znalazłem podejście oparte na zasadach, które moim zdaniem ma dużą zaletę, ponieważ jest koncepcyjnie bardzo jasne i łatwe do dostosowania.

W tym wątku: Jak efektywnie generować losowe macierze korelacji dodatnich i półkolistych? - Opisałem i podałem kod dwóch wydajnych algorytmów generowania macierzy losowej korelacji. Oba pochodzą z artykułu Lewandowskiego, Kurowickiej i Joe (2009), do którego @ssdecontrol przywołał w powyższych komentarzach (wielkie dzięki!).

Proszę zobaczyć moją odpowiedź tam dla wielu liczb, wyjaśnień i kodu Matlab. Tak zwana metoda „winorośli” pozwala na generowanie losowych macierzy korelacji z dowolnym rozkładem korelacji cząstkowych i może być stosowana do generowania macierzy korelacji o dużych wartościach nie przekątnych. Oto przykładowa postać z tego wątku:

Jedyną rzeczą, która zmienia się między wykresami podrzędnymi, jest jeden parametr, który kontroluje, o ile rozkład korelacji cząstkowych jest skoncentrowany wokół .$\pm 1$

Kopiuję mój kod, aby również wygenerować te macierze, aby pokazać, że nie jest on dłuższy niż inne sugerowane tutaj metody. Proszę zobaczyć moją linkowaną odpowiedź dla niektórych wyjaśnień. Wartości betaparampowyższego rysunku wynosiły (a wymiarowość wynosiła 100 ).${50,20,10,5,2,1}$d$100$

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

Aktualizacja: wartości własne

@psarka pyta o wartości własne tych matryc. Na poniższym rysunku wykreślam widma wartości własnych tych samych sześciu macierzy korelacji jak powyżej. Zauważ, że zmniejszają się stopniowo; przeciwnie, metoda sugerowana przez @psarka generalnie daje w wyniku macierz korelacji z jedną dużą wartością własną, ale reszta jest dość jednolita.

Aktualizacja. Naprawdę prosta metoda: kilka czynników

$k<n$$\mathbf W$$k \times n$$\mathbf W \mathbf W^\top$$\mathbf D$$\mathbf B = \mathbf W \mathbf W^\top + \mathbf D$$k={100, 50, 20, 10, 5, 1}$

$k$

Oto kod:

d = 100;    %// number of dimensions
k = 5;      %// number of factors

W = randn(d,k);
S = W*W' + diag(rand(1,d));
S = diag(1./sqrt(diag(S))) * S * diag(1./sqrt(diag(S)));

$a$

import numpy as np
from random import choice
import matplotlib.pyplot as plt

n = 100
a = 2

A = np.matrix([np.random.randn(n) + np.random.randn(1)*a for i in range(n)])
A = A*np.transpose(A)
D_half = np.diag(np.diag(A)**(-0.5))
C = D_half*A*D_half

vals = list(np.array(C.ravel())[0])
plt.hist(vals, range=(-1,1))
plt.show()
plt.imshow(C, interpolation=None)
plt.show()

Hmm, po zrobieniu przykładu w moim języku MatMate widzę, że istnieje już odpowiedź python, co może być preferowane, ponieważ python jest powszechnie używany. Ale ponieważ wciąż masz pytania, pokazuję ci moje podejście przy użyciu języka Matmate-Matrix, być może jest to bardziej samokomunikujące.

Metoda 1
(przy użyciu MatMate):

v=12         // 12 variables
f=3          // subset-correlation based on 3 common factors
vg = v / f   // variables per subsets

 // generate hidden factor-matrix
             // randomu(rows,cols ,lowbound, ubound) gives uniform random matrix 
             //    without explicite bounds the default is: randomu(rows,cols,0,100)
L = {   randomu(vg,f)     || randomu(vg,f)/100  || randomu(vg,f)/100 , _
        randomu(vg,f)/100 || randomu(vg,f)      || randomu(vg,f)/100 , _
        randomu(vg,f)/100 || randomu(vg,f)/100  || randomu(vg,f)     }

 // make sure there is itemspecific variance
 // by appending a diagonal-matrix with random positive entries
L = L || mkdiag(randomu(v,1,10,20)) 
  // make covariance and correlation matrix
cov = L *'   // L multiplied  with its transpose
cor = covtocorr(cov)
                   set ccdezweite=3 ccfeldweite=8
                   list cor
cor = 
   1.000,   0.321,   0.919,   0.489,   0.025,   0.019,   0.019,   0.030,   0.025,   0.017,   0.014,   0.014
   0.321,   1.000,   0.540,   0.923,   0.016,   0.015,   0.012,   0.030,   0.033,   0.016,   0.012,   0.015
   0.919,   0.540,   1.000,   0.679,   0.018,   0.014,   0.012,   0.029,   0.028,   0.014,   0.012,   0.012
   0.489,   0.923,   0.679,   1.000,   0.025,   0.022,   0.020,   0.040,   0.031,   0.014,   0.011,   0.014
   0.025,   0.016,   0.018,   0.025,   1.000,   0.815,   0.909,   0.758,   0.038,   0.012,   0.018,   0.014
   0.019,   0.015,   0.014,   0.022,   0.815,   1.000,   0.943,   0.884,   0.035,   0.012,   0.014,   0.012
   0.019,   0.012,   0.012,   0.020,   0.909,   0.943,   1.000,   0.831,   0.036,   0.013,   0.015,   0.010
   0.030,   0.030,   0.029,   0.040,   0.758,   0.884,   0.831,   1.000,   0.041,   0.017,   0.022,   0.020
   0.025,   0.033,   0.028,   0.031,   0.038,   0.035,   0.036,   0.041,   1.000,   0.831,   0.868,   0.780
   0.017,   0.016,   0.014,   0.014,   0.012,   0.012,   0.013,   0.017,   0.831,   1.000,   0.876,   0.848
   0.014,   0.012,   0.012,   0.011,   0.018,   0.014,   0.015,   0.022,   0.868,   0.876,   1.000,   0.904
   0.014,   0.015,   0.012,   0.014,   0.014,   0.012,   0.010,   0.020,   0.780,   0.848,   0.904,   1.000

Problemem może być tutaj to, że definiujemy bloki podmacierzy, które mają wysokie korelacje wewnątrz, z małą korelacją między nimi, i to nie jest programowo, ale przez ciągłe wyrażenia konkatenacji. Być może to podejście można by modelować bardziej elegancko w Pythonie.

{L = matrix(8,8);  \\ generate an empty factor-loadings-matrix
for(r=1,8, 
   rv=1.0;    \\ remaining variance for variable is 1.0
   for(c=1,8,
        pv=if(c<8,random(100)/100.0,1.0); \\ define randomly part of remaining variance
        cv= pv * rv;  \\ compute current partial variance
        rv = rv - cv;     \\ compute the now remaining variance
        sg = (-1)^(random(100) % 2) ;  \\ also introduce randomly +- signs
        L[r,c] = sg*sqrt(cv) ;  \\ compute factor loading as signed sqrt of cv
       )
     );}

cor = L * L~

a wytworzoną macierzą korelacji jest

     1.000  -0.7111  -0.08648   -0.7806   0.8394  -0.7674   0.6812    0.2765
   -0.7111    1.000   0.06073    0.7485  -0.7550   0.8052  -0.8273   0.05863
  -0.08648  0.06073     1.000    0.5146  -0.1614   0.1459  -0.4760  -0.01800
   -0.7806   0.7485    0.5146     1.000  -0.8274   0.7644  -0.9373  -0.06388
    0.8394  -0.7550   -0.1614   -0.8274    1.000  -0.5823   0.8065   -0.1929
   -0.7674   0.8052    0.1459    0.7644  -0.5823    1.000  -0.7261   -0.4822
    0.6812  -0.8273   -0.4760   -0.9373   0.8065  -0.7261    1.000   -0.1526
    0.2765  0.05863  -0.01800  -0.06388  -0.1929  -0.4822  -0.1526     1.000

Być może generuje to macierz korelacji z dominującymi składnikami głównymi z powodu skumulowanej reguły generowania dla macierzy ładunków czynnikowych. Może być również lepiej zapewnić pozytywną definitywność, czyniąc ostatnią część wariancji unikalnym czynnikiem. Zostawiłem to w programie, aby skupić się na ogólnej zasadzie.

Macierz korelacji 100 x 100 miała następujące częstotliwości korelacji (w zaokrągleniu do 1 dec)

    e    f            e: entry(rounded) f: frequency
  -----------------------------------------------------
  -1.000, 108.000
  -0.900, 460.000
  -0.800, 582.000
  -0.700, 604.000
  -0.600, 548.000
  -0.500, 540.000
  -0.400, 506.000
  -0.300, 482.000
  -0.200, 488.000
  -0.100, 464.000
   0.000, 434.000
   0.100, 486.000
   0.200, 454.000
   0.300, 468.000
   0.400, 462.000
   0.500, 618.000
   0.600, 556.000
   0.700, 586.000
   0.800, 536.000
   0.900, 420.000
   1.000, 198.000

[aktualizacja]. Hmm, matryca 100 x 100 jest źle kondycjonowana; Pari / GP nie może poprawnie określić wartości własnych za pomocą polroots (charpoly ()) - funkcja nawet z dokładnością do 200 cyfr. Zrobiłem obrót Jacobiego, aby utworzyć formę pca na macierzy obciążeń L i znajdowałem głównie bardzo małe wartości własne, wydrukowałem je w logarytmach do podstawy 10 (które podają z grubsza pozycję punktu dziesiętnego). Czytaj od lewej do prawej, a następnie wiersz po rzędzie:

log_10(eigenvalues):
   1.684,   1.444,   1.029,   0.818,   0.455,   0.241,   0.117,  -0.423,  -0.664,  -1.040
  -1.647,  -1.799,  -1.959,  -2.298,  -2.729,  -3.059,  -3.497,  -3.833,  -4.014,  -4.467
  -4.992,  -5.396,  -5.511,  -6.366,  -6.615,  -6.834,  -7.535,  -8.138,  -8.263,  -8.766
  -9.082,  -9.482,  -9.940, -10.167, -10.566, -11.110, -11.434, -11.788, -12.079, -12.722
 -13.122, -13.322, -13.444, -13.933, -14.390, -14.614, -15.070, -15.334, -15.904, -16.278
 -16.396, -16.708, -17.022, -17.746, -18.090, -18.358, -18.617, -18.903, -19.186, -19.476
 -19.661, -19.764, -20.342, -20.648, -20.805, -20.922, -21.394, -21.740, -21.991, -22.291
 -22.792, -23.184, -23.680, -24.100, -24.222, -24.631, -24.979, -25.161, -25.282, -26.211
 -27.181, -27.626, -27.861, -28.054, -28.266, -28.369, -29.074, -29.329, -29.539, -29.689
 -30.216, -30.784, -31.269, -31.760, -32.218, -32.446, -32.785, -33.003, -33.448, -34.318

[aktualizacja 2]
Metoda 2 (b)
Ulepszenie może polegać na zwiększeniu wariancji specyficznej dla przedmiotu do pewnego nie-marginalnego poziomu i zredukowaniu do rozsądnie mniejszej liczby wspólnych czynników (na przykład liczba całkowita-kwadratowa numeru pozycji):

{  dimr = 100;
   dimc = sqrtint(dimr);        \\ 10 common factors
   L = matrix(dimr,dimr+dimc);  \\ loadings matrix 
                                \\     with dimr itemspecific and 
                                \\          dimc common factors
   for(r=1,dim, 
         vr=1.0;                \\ complete variance per item 
         vu=0.05+random(100)/1000.0;   \\ random variance +0.05
                                       \\ for itemspecific variance
         L[r,r]=sqrt(vu);              \\ itemspecific factor loading  
         vr=vr-vu;
         for(c=1,dimc,
                cv=if(c<dimc,random(100)/100,1.0)*vr;
                vr=vr-cv;
                L[r,dimr+c]=(-1)^(random(100) % 2)*sqrt(cv)
             )
        );}

   cov=L*L~
   cp=charpoly(cov)   \\ does not work even with 200 digits precision
   pr=polroots(cp)    \\ spurious negative and complex eigenvalues...

Struktura wyniku

pod względem rozkładu korelacji:

pozostaje podobny (również nieprzyjemny brak rozkładu przez PariGP), ale wartości własne, gdy zostały znalezione przez rotację jacobi macierzy loadings, mają teraz lepszą strukturę, dla nowo obliczonego przykładu otrzymałem wartości własne jako

log_10(eigenvalues):
   1.677,   1.326,   1.063,   0.754,   0.415,   0.116,  -0.262,  -0.516,  -0.587,  -0.783
  -0.835,  -0.844,  -0.851,  -0.854,  -0.858,  -0.862,  -0.862,  -0.868,  -0.872,  -0.873
  -0.878,  -0.882,  -0.884,  -0.890,  -0.895,  -0.896,  -0.896,  -0.898,  -0.902,  -0.904
  -0.904,  -0.909,  -0.911,  -0.914,  -0.920,  -0.923,  -0.925,  -0.927,  -0.931,  -0.935
  -0.939,  -0.939,  -0.943,  -0.948,  -0.951,  -0.955,  -0.956,  -0.960,  -0.967,  -0.969
  -0.973,  -0.981,  -0.986,  -0.989,  -0.997,  -1.003,  -1.005,  -1.011,  -1.014,  -1.019
  -1.022,  -1.024,  -1.031,  -1.038,  -1.040,  -1.048,  -1.051,  -1.061,  -1.064,  -1.068
  -1.070,  -1.074,  -1.092,  -1.092,  -1.108,  -1.113,  -1.120,  -1.134,  -1.139,  -1.147
  -1.150,  -1.155,  -1.158,  -1.166,  -1.171,  -1.175,  -1.184,  -1.184,  -1.192,  -1.196
  -1.200,  -1.220,  -1.237,  -1.245,  -1.252,  -1.262,  -1.269,  -1.282,  -1.287,  -1.290

$A$$B$$\lambda A + (1-\lambda)B$$\lambda$

$A$$B$$C$$\lambda_A A + \lambda_B B + \lambda_C C$$\sum \lambda = 1$$\lambda \geq 0$

R ma pakiet (klasterGeneration), który implementuje metodę w:

• Joe, H. (2006) Generowanie macierzy losowej korelacji na podstawie korelacji częściowych . Journal of Multivariate Analysis, 97, 2177--2189.

Przykład:

> (cormat10 = clusterGeneration::rcorrmatrix(10, alphad = 1/100000000000000))
        [,1]   [,2]    [,3]     [,4]     [,5]   [,6]   [,7]    [,8]     [,9]   [,10]
 [1,]  1.000  0.344 -0.1406 -0.65786 -0.19411  0.246  0.688 -0.6146  0.36971 -0.1052
 [2,]  0.344  1.000 -0.4256 -0.35512  0.15973  0.192  0.340 -0.4907 -0.30539 -0.6104
 [3,] -0.141 -0.426  1.0000  0.01775 -0.61507 -0.485 -0.273  0.3492 -0.30284  0.1647
 [4,] -0.658 -0.355  0.0178  1.00000  0.00528 -0.335 -0.124  0.5256 -0.00583 -0.0737
 [5,] -0.194  0.160 -0.6151  0.00528  1.00000  0.273 -0.350 -0.0785  0.08285  0.0985
 [6,]  0.246  0.192 -0.4847 -0.33531  0.27342  1.000  0.278 -0.2220 -0.11010  0.0720
 [7,]  0.688  0.340 -0.2734 -0.12363 -0.34972  0.278  1.000 -0.6409  0.40314 -0.2800
 [8,] -0.615 -0.491  0.3492  0.52557 -0.07852 -0.222 -0.641  1.0000 -0.50796  0.1461
 [9,]  0.370 -0.305 -0.3028 -0.00583  0.08285 -0.110  0.403 -0.5080  1.00000  0.3219
[10,] -0.105 -0.610  0.1647 -0.07373  0.09847  0.072 -0.280  0.1461  0.32185  1.0000
> cormat10[lower.tri(cormat10)] %>% psych::describe()
   vars  n  mean   sd median trimmed mad   min  max range skew kurtosis   se
X1    1 45 -0.07 0.35  -0.08   -0.07 0.4 -0.66 0.69  1.35 0.03       -1 0.05

Niestety nie wydaje się możliwe symulowanie korelacji, które następują po tym rozkładzie równomiernym. Wydaje się, że silniej koreluje, kiedyalphad są ustawione na bardzo małe wartości, ale nawet przy 1/100000000000000tym zakres korelacji wzrósłby tylko do około 1,40.

Niemniej jednak mam nadzieję, że może to komuś się przydać.