Wie können zufällige positiv-semidefinite Korrelationsmatrizen effizient generiert werden?

38

Ich möchte in der Lage sein, Korrelationsmatrizen mit positivem Semidefinit (PSD) effizient zu erzeugen. Meine Methode verlangsamt sich dramatisch, wenn ich die zu generierenden Matrizen vergrößere.

  1. Können Sie effiziente Lösungen vorschlagen? Wenn Sie Beispiele in Matlab kennen, wäre ich Ihnen sehr dankbar.
  2. Wie würden Sie beim Generieren einer PSD-Korrelationsmatrix die Parameter auswählen, um die zu generierenden Matrizen zu beschreiben? Eine durchschnittliche Korrelation, Standardabweichung von Korrelationen, Eigenwerte?
Eduardas
quelle

Antworten:

16

Sie können es rückwärts machen: Jede Matrix (die Menge aller symmetrischen PSD-Matrizen) kann als zerlegt werden p × pCR++pp×p

OC=OTDO wobei eine orthonormale Matrix istO

Um , generieren Sie zuerst eine Zufallsbasis (wobei Zufallsvektoren sind, typischerweise in ). Verwenden Sie von dort den Gram-Schmidt-Orthogonalisierungsprozess, um( v 1 , . . . , V p ) v i ( - 1 , 1 ) ( u 1 , . . . . , U p ) = OO(v1,...,vp)vi(1,1)(u1,....,up)=O

R verfügt über eine Reihe von Paketen, mit denen die GS-Orthogonalisierung auf zufälliger Basis effizient durchgeführt werden kann, auch für große Dimensionen, z. B. das "ferne" Paket. Obwohl Sie den GS-Algorithmus im Wiki finden, ist es wahrscheinlich besser, das Rad nicht neu zu erfinden und sich für eine Matlab-Implementierung zu entscheiden (es gibt sicherlich eine, ich kann keine empfehlen).

Schließlich ist eine Diagonalmatrix, deren Elemente alle positiv sind (dies ist wiederum leicht zu erzeugen: Zufallszahlen erzeugen , quadrieren, sortieren und in die Diagonale einer Identität durch Matrix platzieren).p p pDppp

user603
quelle
3
(1) Beachten Sie, dass das resultierende keine Korrelationsmatrix ist (wie vom OP gefordert), da es keine auf der Diagonale hat. Natürlich kann die Skalierung geändert werden, indem die Diagonale auf , wobei eine Diagonalmatrix mit derselben Diagonale wie . (2) Wenn ich mich nicht irre, wird dies zu Korrelationsmatrizen führen, bei denen alle nicht diagonalen Elemente um konzentriert sind , so dass es keine Flexibilität gibt, nach der das OP gesucht hat (OP wollte in der Lage sein, eine durchschnittliche Korrelation zu setzen) , Standardabweichung von Korrelationen, Eigenwerten " )E - 1 / 2 C E - 1 / 2 E C 0CE1/2CE1/2EC0
Amöbe sagt Reinstate Monica
@amoeba: Ich werde auf (2) eingehen, da die Lösung für (1), wie Sie hervorheben, trivial ist. Die Einzahlcharakterisierung der 'Form' (die Beziehung zwischen den In- und Out-Diagonalelementen) einer PSD-Matrix (und damit einer Kovarianz- und auch einer Korrelationsmatrix) ist ihre Bedingungszahl. Und das oben beschriebene Verfahren ermöglicht die vollständige Kontrolle darüber. Die "Konzentration der Elemente außerhalb der Diagonale um 0" ist kein Merkmal der Methode zur Erzeugung von PSD-Matrizen, sondern eine Folge der Einschränkungen, die erforderlich sind, um sicherzustellen, dass die Matrix PSD ist, und der Tatsache, dass groß ist. p
User603
Wollen Sie damit sagen, dass alle großen PSD-Matrizen nicht diagonale Elemente nahe Null haben? Ich stimme nicht zu, es ist nicht so. Überprüfen Sie hier meine Antwort auf einige Beispiele: Wie man eine zufällige Korrelationsmatrix erzeugt, die ungefähr normalverteilte Einträge außerhalb der Diagonale mit gegebener Standardabweichung aufweist? Aber man kann direkt sehen, dass dies nicht der Fall ist, weil eine quadratische Matrix mit allen auf der Diagonale und einem festen Wert überall außerhalb der Diagonale PSD ist und beliebig groß sein kann (aber natürlich unter ). ρ 1ρρ1
Amöbe sagt Reinstate Monica
@amoeba: Dann habe ich zu Unrecht angenommen, dass die Off-Diagonale großer Korrelationsmatrizen, wenn sie sowohl positiv als auch negativ sein dürfen, nahe bei 0 liegt. Danke für das aufschlussreiche Beispiel.
User603
1
Ich habe eine sehr schöne Abhandlung über das Erzeugen zufälliger Korrelationsmatrizen gelesen und hier meine eigene Antwort gegeben (sowie eine andere Antwort in diesem verknüpften Thread). Ich denke, du findest es vielleicht interessant.
Amöbe sagt Reinstate Monica
27

Ein Artikel zur Generierung zufälliger Korrelationsmatrizen auf der Basis von Reben und einer Methode mit erweiterten Zwiebeln von Lewandowski, Kurowicka und Joe (LKJ), 2009, bietet eine einheitliche Behandlung und Darstellung der beiden effizienten Methoden zur Generierung zufälliger Korrelationsmatrizen. Beide Methoden ermöglichen die Erzeugung von Matrizen aus einer gleichmäßigen Verteilung in einem bestimmten, nachstehend definierten Sinne, sind einfach zu implementieren, schnell und haben den zusätzlichen Vorteil, dass sie amüsante Namen haben.

Eine reelle symmetrische Matrix mit der Größe mit Einsen auf der Diagonale hat eindeutige Elemente außerhalb der Diagonale und kann daher als Punkt in parametrisiert werden. . Jeder Punkt in diesem Raum entspricht einer symmetrischen Matrix, aber nicht alle von ihnen sind positiv-definit (wie Korrelationsmatrizen sein müssen). Korrelationsmatrizen bilden daher eine Teilmenge von (eigentlich eine zusammenhängende konvexe Teilmenge), und beide Methoden können Punkte aus einer gleichmäßigen Verteilung über diese Teilmenge erzeugen.d ( d - 1 ) / 2 R d ( d - 1 ) / 2 R d ( d - 1 ) / 2d×dd(d1)/2Rd(d1)/2Rd(d1)/2

Ich werde meine eigene MATLAB-Implementierung jeder Methode bereitstellen und sie mit veranschaulichen .d=100


Zwiebelmethode

Die Zwiebelmethode stammt aus einer anderen Veröffentlichung (Lit. 3 in LKJ) und hat ihren Namen der Tatsache zu verdanken, dass die Korrelationsmatrizen beginnend mit Matrix und spaltenweise und zeilenweise wachsen. Die resultierende Verteilung ist gleichmäßig. Ich verstehe die Mathematik hinter der Methode nicht wirklich (und bevorzuge sowieso die zweite Methode), aber hier ist das Ergebnis:1×1

Zwiebelmethode

Hier und unter dem Titel jedes Unterplots werden der kleinste und der größte Eigenwert sowie die Determinante (Produkt aller Eigenwerte) angezeigt. Hier ist der Code:

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

Erweiterte Zwiebelmethode

LKJ modifizieren diese Methode geringfügig, um Korrelationsmatrizen aus einer Verteilung abtasten zu können, die proportional zu [ d e t istC[detC]η1ηη=1η=1,10,100,1000,10000,100000

Erweiterte Zwiebelmethode

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

Rebe-Methode

d(d1)/2[1,1]η[detC]η1

Rebe-Methode

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

Weinbaumethode mit manueller Probenahme von Teilkorrelationen

±1[0,1][1,1]α=β=50,20,10,5,2,1. Je kleiner die Parameter der Beta-Verteilung sind, desto stärker ist sie in Randnähe konzentriert.

Rebmethode mit manueller Probenahme

Beachten Sie, dass in diesem Fall die Verteilung nicht garantiert permutationsinvariant ist, sodass ich Zeilen und Spalten nach der Generierung zusätzlich zufällig permutiere.

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

So sehen die Histogramme der nicht diagonalen Elemente für die obigen Matrizen aus (Varianz der Verteilung steigt monoton an):

Nicht diagonale Elemente


Update: unter Verwendung von Zufallsfaktoren

k<dWk×dWWDB=WW+DC=E1/2BE1/2EBk=100,50,20,10,5,1

zufällige Korrelationsmatrizen aus zufälligen Faktoren

Und der Code:

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

Hier ist der Umhüllungscode, der zum Generieren der Zahlen verwendet wird:

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
Amöbe sagt Reinstate Monica
quelle
2
Das ist ein fantastischer Abriss, ich bin froh, dass ich etwas gesagt habe!
Shadowtalker
Als ich den Matlab-Code für die rebenbasierte Korrelationsmatrix in R übersetzt und getestet habe, war die Korrelationsdichte in Spalte 1 immer anders als in späteren Spalten. Es kann sein, dass ich etwas falsch übersetzt habe, aber vielleicht hilft dieser Hinweis jemandem.
Charlie
3
Für R-Benutzer implementiert die Funktion rcorrmatrix im Paket clusterGeneration (geschrieben von W Qui und H. Joe) die vine-Methode.
RNM
15

AATAyT(ATA)y0yyT(ATA)y=(Ay)TAy=||Ay||Das ist nicht negativ. Versuchen Sie es in Matlab einfach

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

Je nach Anwendung erhalten Sie möglicherweise nicht die gewünschte Verteilung der Eigenwerte. Kwaks Antwort ist in dieser Hinsicht viel besser. Die Eigenwerte Xdieses Code-Snippets sollten der Marchenko-Pastur-Verteilung folgen.

Um die Korrelationsmatrizen von Aktien zu simulieren, möchten Sie vielleicht einen etwas anderen Ansatz:

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
quelle
1
Wären Sie bereit, Ihre Funktion ascii_hist zufällig mit anderen zu teilen?
Stadt
@btown der Rand ist zu klein, um es zu enthalten!
Shabbychef
1
yT(ATA)y=(Ay)TAy=||Ay||
8

DA=QDQTQ

JM ist kein Statistiker
quelle
M.: Netter Hinweis: Dies scheint die effizienteste Lösung zu sein (asymptotisch).
Whuber
3
@whuber: Heh, ich habe es von Golub und Van Loan (natürlich) abgeholt; Ich nutze dies die ganze Zeit, um Testmatrizen für das Stresstesten von Eigenwert- / Singularwert-Routinen zu generieren. Wie aus dem Artikel hervorgeht, entspricht dies im Wesentlichen der QR-Zerlegung einer Zufallsmatrix, wie sie von Kwak vorgeschlagen wurde, mit der Ausnahme, dass sie effizienter durchgeführt wird. Es gibt eine MATLAB-Implementierung davon in der Text Matrix Toolbox von Higham, BTW.
JM ist kein Statistiker
M .:> Danke für die Matlab-Implementierung. Würdest du zufällig einen Haar-Pseudozufallsmatrixgenerator in R kennen?
user603
@kwak: Keine Ahnung, aber wenn es noch keine Implementierung gibt, sollte es nicht allzu schwierig sein, den MATLAB-Code in R zu übersetzen (ich kann versuchen, einen hochzuspielen, wenn es wirklich keinen gibt). die einzige voraussetzung ist ein anständiger generator für pseudozufällige normale variaten, die ich mir sicher habe r.
JM ist kein Statistiker
M .:> Ja, ich werde es wahrscheinlich selbst übersetzen. Danke für die Links, Best.
user603
4

Sie haben keine Verteilung für die Matrizen angegeben. Zwei gebräuchliche sind die Wishart- und die inverse Wishart-Verteilung. Die Bartlett-Zerlegung ergibt eine Cholesky-Faktorisierung einer zufälligen Wishart-Matrix (die auch effizient gelöst werden kann, um eine zufällige inverse Wishart-Matrix zu erhalten).

Tatsächlich ist der Cholesky-Raum eine bequeme Möglichkeit, andere Arten von zufälligen PSD-Matrizen zu generieren, da Sie nur sicherstellen müssen, dass die Diagonale nicht negativ ist.

Simon Byrne
quelle
> Nicht zufällig: Zwei von demselben Whishard erzeugte Matrizen sind nicht unabhängig voneinander. Wenn Sie vorhaben, den Whishart bei jeder Generation zu ändern, wie beabsichtigen Sie dann, diesen Whishart überhaupt zu erzeugen?
user603
@kwak: Ich verstehe Ihre Frage nicht: Die Bartlett-Zerlegung ergibt unabhängige Draws aus derselben Wishart-Distribution.
Simon Byrne
> Lassen Sie mich das umformulieren. Woher erhalten Sie die Skalierungsmatrix Ihrer Whishart-Verteilung?
user603
1
@kwak: Es ist ein Parameter der Distribution und daher fest. Sie wählen es zu Beginn basierend auf den gewünschten Eigenschaften Ihrer Distribution (wie dem Mittelwert) aus.
Simon Byrne
3

UTSU

gappy
quelle
Wenn die Einträge aus einer Normalverteilung und nicht aus einer Uniform generiert werden, sollte die von Ihnen erwähnte Zerlegung SO (n) -invariant sein (und daher im Verhältnis zum Haar-Maß gleich verteilt sein).
whuber
Interessant. Können Sie hierfür eine Referenz angeben?
gappy
1
> das problem bei dieser methode ist, dass sie das verhältnis von kleinstem zu größtem eigenwert nicht steuern können (und ich denke, dass, wenn die größe ihres zufällig generierten datensatzes unendlich wird, dieses verhältnis gegen 1 konvergiert).
user603
1

Wenn Sie mehr Kontrolle über Ihre generierte symmetrische PSD-Matrix haben möchten, z. B. einen synthetischen Validierungsdatensatz generieren möchten, stehen Ihnen eine Reihe von Parametern zur Verfügung. Eine symmetrische PSD-Matrix entspricht einer Hyperellipse im N-dimensionalen Raum mit allen damit verbundenen Freiheitsgraden:

  1. Drehungen.
  2. Längen der Achsen.

Für eine zweidimensionale Matrix (dh eine 2D-Ellipse) haben Sie also 1 Drehung + 2 Achsen = 3 Parameter.

Σ=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
quelle
0

Ein billiger und fröhlicher Ansatz, den ich zum Testen verwendet habe, besteht darin, m N (0,1) n-Vektoren V [k] zu generieren und dann P = d * I + Summe {V [k] * V [k] '} zu verwenden. als nxn psd matrix. Mit m <n ist dies für d = 0 singulär und für kleines d hat es eine hohe Bedingungszahl.


quelle
2
> das problem bei dieser methode ist, dass sie das verhältnis von kleinstem zu größtem eigenwert nicht steuern können (und ich denke, dass, wenn die größe ihres zufällig generierten datensatzes unendlich wird, dieses verhältnis gegen 1 konvergiert).
user603
> Außerdem ist das Verfahren nicht sehr effizient (von computationnal Sicht)
user603
1
Ihre "Zufallsmatrix" ist eine speziell strukturierte Matrix mit der Bezeichnung "Diagonale plus Rang-1-Matrix" (DR1-Matrix), also keine wirklich gute repräsentative Zufallsmatrix.
JM ist kein Statistiker