Wie werden korrelierte Zufallszahlen generiert (gegebene Mittelwerte, Varianzen und Grad der Korrelation)?

53

Es tut mir leid, wenn dies ein bisschen zu grundlegend erscheint, aber ich schätze, ich versuche hier nur, das Verständnis zu bestätigen. Ich habe das Gefühl, dass ich dies in zwei Schritten tun müsste, und ich habe angefangen, Korrelationsmatrizen zu erstellen, aber es scheint erst sehr involviert zu sein. Ich suche eine prägnante Erklärung (idealerweise mit Hinweisen auf eine Pseudocodelösung) für einen guten, idealerweise schnellen Weg, um korrelierte Zufallszahlen zu generieren.

Angesichts zweier Pseudozufallsvariablen Größe und Gewicht mit bekannten Mitteln und Varianzen sowie einer gegebenen Korrelation versuche ich im Grunde zu verstehen, wie dieser zweite Schritt aussehen sollte:

   height = gaussianPdf(height.mean, height.variance)
   weight = gaussianPdf(correlated_mean(height.mean, correlation_coefficient), 
                        correlated_variance(height.variance, 
                        correlation_coefficient))
  • Wie berechne ich den korrelierten Mittelwert und die Varianz? Aber ich möchte bestätigen, dass das hier wirklich das relevante Problem ist.
  • Muss ich auf Matrixmanipulation zurückgreifen? Oder habe ich noch etwas sehr Falsches in meiner grundlegenden Herangehensweise an dieses Problem?
Joseph Weissman
quelle
1
Ich bin mir nicht sicher, ob ich Sie richtig verstehe, aber Sie müssen nicht den "korrelierten Mittelwert und die Varianz" berechnen. Wenn Sie davon ausgehen, dass die Variablen bivariat normal sind, sollten Sie die einzelnen Mittelwerte und Varianzen sowie die Korrelation angeben. Gibt es eine bestimmte Software, die Sie dafür verwenden möchten?
mark999
1
Außerdem: Wie kann ich Daten mit einer vorgegebenen Korrelationsmatrix generieren?
gung - Wiedereinsetzung von Monica

Antworten:

44

Zur Beantwortung Ihrer Frage "Ein guter, idealerweise schneller Weg, um korrelierte Zufallszahlen zu generieren": Bei einer gewünschten Varianz-Kovarianz-Matrix , die per Definition eindeutig positiv ist, ist ihre Cholesky-Zerlegung: = ; ist die untere Dreiecksmatrix.C L L T LCCLLTL

Wenn Sie nun mit dieser Matrix einen unkorrelierten Zufallsvariablenvektor projizieren, ist die resultierende Projektion die von korrelierten Zufallsvariablen.X Y = L XLXY=LX

Eine kurze Erklärung dafür finden Sie hier .

usεr11852 sagt Reinstate Monic
quelle
Vielen Dank! Das war enorm hilfreich. Ich glaube, ich habe zumindest ein besseres Gespür dafür, was ich als Nächstes betrachten muss.
Joseph Weissman
7
Gilt diese Methode nur für Gauß-Verteilungen (wie in der Frage angegeben) oder kann sie zum Generieren von korrelierten Variablen verwendet werden, die anderen Verteilungen folgen? Wenn nicht, kennen Sie eine Methode, die in diesem Fall verwendet werden könnte?
user000001
1
@ Michael: Ja. Vorausgesetzt, dass eine gültige Kovarianzmatrix ist, ist die Cholesky-Zerlegung der schnellste Weg. Sie könnten auch die (symmetrische) Quadratwurzel- Matrix von mit SVD erhalten (also , wobei von ), aber das wäre mehr auch teuer. X C C = X X = X X T X = U S 0.5 V T C = U S V TCXCC=XX=XXTX=US0.5VTC=USVT
usεr11852 sagt Reinstate Monic
1
@ Michael: Natürlich. Ihre Kovarianz wird (ungefähr) gleich sein, nicht die Zahlen selbst.
usεr11852 sagt Reinstate Monic
1
@Sid: Jede kontinuierliche Verteilung, die nicht auf der gesamten realen Leitung unterstützt wird, schlägt sofort fehl. Wenn wir zum Beispiel ein einheitliches , können wir nicht garantieren, dass die "korrelierten Zahlen" in ; Ähnliches gilt für ein Poisson, wenn wir nicht diskrete Zahlen erhalten. Darüber hinaus schlägt auch jede Verteilung fehl , bei der die Summe der Verteilungen noch nicht dieselbe Verteilung ist (z. B. die Summierung der Verteilung führt nicht zu Verteilungen ). In allen genannten Fällen werden die produzierten Zahlen mit dem korreliert , sie entsprechen jedoch nicht der von uns gestarteten Verteilung. [ 0 , 1 ] t t CU[0,1][0,1]ttC
usεr11852 sagt Reinstate Monic
36

+1 an @ user11852 und @ jem77bfp, das sind gute Antworten. Lassen Sie mich dies aus einer anderen Perspektive betrachten, nicht weil ich denke, dass es in der Praxis unbedingt besser ist , sondern weil ich es für lehrreich halte. Hier sind einige relevante Fakten, die wir bereits kennen:

  1. X Y N ( 0 , 1 )r ist die Steigung der Regressionslinie , wenn beide und sind standardisiert , dh , XYN(0,1)
  2. Y Xr2 ist der Anteil der Varianz in der auf die Varianz in , YX



    (auch aus den Regeln für Abweichungen ):

  3. Die Varianz einer Zufallsvariablen multipliziert mit einer Konstanten ist das Quadrat der ursprünglichen Varianz:
    Var[aX]=a2Var[X]
  4. Varianzen addieren , dh die Varianz der Summe zweier Zufallsvariablen (sofern sie unabhängig sind) ist die Summe der beiden Varianzen:
    Var[X+ε]=Var[X]+Var[ε]

Jetzt können wir diese vier Fakten kombinieren, um zwei normale Standardvariablen zu erstellen, deren Populationen eine bestimmte Korrelation haben, (genauer gesagt ), obwohl die von Ihnen generierten Stichproben unterschiedliche Stichproben-Korrelationen haben. Die Idee ist, eine Pseudozufallsvariable zu erstellen , die normal ist, , und dann einen Koeffizienten und eine Fehlervarianz , so dass , wobei . (Beachten Sie, dass muss, damit dies funktioniert, und dass außerdem .) Sie beginnen also mit demρ X N ( 0 , 1 ) a v e Y ~ N ( 0 , a 2 + v e ) a 2 + v e = 1 | a | 1 a = r r a 1 - r 2 x i e i v e y irρXN(0,1)aveYN(0,a2+ve)a2+ve=1|a| 1a=rr das du willst; das ist dein Koeffizient, . Dann ermitteln Sie die benötigte Fehlervarianz: . (Wenn Ihre Software die Verwendung der Standardabweichung erfordert, nehmen Sie die Quadratwurzel dieses Werts.) Generieren Sie schließlich für jede von Ihnen generierte Pseudozufallsvariable eine Pseudozufallsvariable mit der entsprechenden Fehlervarianz . und die korrelierte Pseudozufallsvariable durch Multiplizieren und Addieren zu berechnen . a1r2xieiveyi

Wenn Sie dies in R tun möchten, könnte der folgende Code für Sie funktionieren:

correlatedValue = function(x, r){
  r2 = r**2
  ve = 1-r2
  SD = sqrt(ve)
  e  = rnorm(length(x), mean=0, sd=SD)
  y  = r*x + e
  return(y)
}

set.seed(5)
x = rnorm(10000)
y = correlatedValue(x=x, r=.5)

cor(x,y)
[1] 0.4945964

(Edit: Ich habe vergessen zu erwähnen :) Wie ich es beschrieben habe, liefert diese Prozedur zwei normale korrelierte Standardvariablen. Wenn Sie keine Standardnormalen möchten , aber möchten, dass die Variablen bestimmte Mittelwerte (nicht 0) und SDs (nicht 1) haben, können Sie sie transformieren, ohne die Korrelation zu beeinträchtigen. Sie subtrahieren also den beobachteten Mittelwert, um sicherzustellen, dass der Mittelwert genau ist. Multiplizieren Sie die Variable mit der gewünschten SD und addieren Sie dann den gewünschten Mittelwert. Wenn Sie möchten, dass der beobachtete Mittelwert normal um den gewünschten Mittelwert schwankt, würden Sie die anfängliche Differenz zurückaddieren. Dies ist im Wesentlichen eine umgekehrte Z-Score-Transformation. Da es sich um eine lineare Transformation handelt, hat die transformierte Variable dieselbe Korrelation mit der anderen Variablen wie zuvor. 0

Auch hier können Sie in der einfachsten Form nur ein Paar korrelierter Variablen generieren (dies könnte skaliert werden, wird aber schnell hässlich) und ist sicherlich nicht die bequemste Methode, um die Aufgabe zu erledigen. In R möchten Sie ? Mvrnorm im MASS- Paket verwenden, da dies einfacher ist und Sie mit einer bestimmten Populationskorrelationsmatrix viele Variablen generieren können. Dennoch finde ich es lohnend, diesen Prozess durchlaufen zu haben, um zu sehen, wie sich einige Grundprinzipien auf einfache Weise auswirken.

gung - Wiedereinsetzung von Monica
quelle
Dieser im Wesentlichen regressive Ansatz ist besonders nützlich, wenn man ein zufälliges Y erzeugen lässt, das mit einer beliebigen Anzahl vorhandener X- "Prädiktoren" korreliert . Habe ich in diesem Verständnis recht?
TTNPHNS
Es hängt genau davon ab, welches Korrelationsmuster zwischen den gewünschten Variablen (@ttnphns) vorliegt. Sie können dies nacheinander durchlaufen, aber es würde langweilig werden. Um viele korrelierte Variablen mit einem bestimmten Muster zu erstellen, ist es besser, die Cholesky-Zerlegung zu verwenden.
gung - Wiedereinsetzung von Monica
Wissen Sie, wie Sie mit Cholesky ein Y-Korrelat (ungefähr wie in Ihrer Methode) anhand eines Korrelationsvektors mit mehreren vorhandenen (nicht simulierten) X erzeugen können?
TTNPHNS
@ttnphns, möchten Sie ein einzelnes Y mit einer gegebenen Populationskorrelation mit einer Menge von X erzeugen, nicht mit einer Menge von p Variablen, die alle vorgegebene Populationskorrelationen haben? Eine einfache Möglichkeit wäre, eine Regressionsgleichung zu schreiben, um einen einzelnen Y-Hut aus Ihren X zu generieren, und dann die obige Methode zu verwenden, um Y als Korrelat Ihres Y-Huts zu generieren. Sie können eine neue Frage stellen, wenn Sie möchten.
gung - Wiedereinsetzung von Monica
1
Das habe ich in meinem ersten Kommentar gemeint: Diese Methode ist eine direkte Erweiterung dessen, wovon Sie in Ihrer Antwort sprechen: im Wesentlichen eine Regressionsmethode (Hat-Methode).
TTNPHNS
16

Im Allgemeinen keine einfache Sache zu tun, aber ich glaube , dass es Pakete für sind multivariate Normal variable Generation (zumindest in R finden Sie mvrnormim MASSPaket), in dem gerade eingegebenen eine Kovarianzmatrix und eine mittlere Vektor.

(X1,X2)F(x1,x2)Fx2

FX1(x1)=F(x1,x2)dx2.
FX11FX1ξ1[0,1]x^1=FX11(ξ)

Da wir nun eine Koordinate haben, müssen wir sie in unsere ursprüngliche Verteilungsfunktion und dann eine bedingte Verteilungsfunktion mit der Bedingung : wobei eine Wahrscheinlichkeitsdichtefunktion von ist marginale Verteilung; dh .F(x1,x2)x1=x^1

F(x2|X1=x^1)=F(x^1,x2)fX1(x^1),
fX1X1FX1(x1)=fX1(x1)

Dann erzeugen Sie wieder eine gleichmäßig verteilte Variable auf (unabhängig von ) und fügen sie in die Umkehrung von . Daher erhalten Sie ; das heißt, erfüllt . Diese Methode kann auf Vektoren mit mehr Dimensionen verallgemeinert werden, hat jedoch den Nachteil, dass Sie viele Funktionen analytisch oder numerisch berechnen müssen. Die Idee finden Sie auch in diesem Artikel: http://www.econ-pol.unisi.it/dmq/pdf/DMQ_WP_34.pdf .ξ2[0,1]ξ1F(x2|X1=x^1)x^2=(F(x2|X1=x^1))1(ξ)x^2F(x^2|X1=x^1)=ξ

Wenn Sie die Bedeutung des Einfügens einer einheitlichen Variablen in eine inverse Wahrscheinlichkeitsverteilungsfunktion nicht verstehen, versuchen Sie, eine Skizze des univariaten Falls zu erstellen, und erinnern Sie sich dann an die geometrische Interpretation der inversen Funktion.

jem77bfp
quelle
Kluge Idee! Hat einfache intuitive Anziehungskraft. Aber ja scheint rechenintensiv.
MichaelChirico
(+1) sehr guter Punkt. Am Anfang wäre es besser zu sagen, dass , dann fließt es natürlicher, zuerst eine mariginal Distribution zu generieren und dann die bedingte Verteilung. Sehr gut! fX,Y(x,y)=fX(x)fY|X(y)
KevinKim
1

Wenn Sie bereit sind, auf Effizienz zu verzichten, können Sie einen Wegwerfalogorithmus verwenden. Ihr Vorteil ist, dass sie beliebige Distributionen (nicht nur Gauß) zulässt.

Beginnen Sie, indem Sie zwei unkorrelierte Folgen von Zufallszahlen und mit beliebigen Verteilungen erzeugen . Sei durch den gewünschten Wert des Korrelationskoeffizienten. Dann machen Sie folgendes:{xi}i=1N{yi}i=1NC

1) Berechne den Korrelationskoeffizientencold=corr({xi},{yi})

2) Generiere zwei Zufallszahlen undn1n2:1n1,2N

3) Tausche die Zahlen undxn1xn2

4) Neue Korrelationcnew=corr({xi},{yi})

5) Wenndann den Swap behalten. Andernfalls machen Sie den Swap rückgängig.|Ccnew|<|Ccold|

6) Wenn stop, sonst gehe zu 1)|Cc|<ϵ

Zufällige Swaps verändern die marginale Verteilung von .xi

Viel Glück!

F. Jatpil
quelle
xicorr(xi,yi)
xi{xi}ycorr(xi,yi)corr({xi},{yi})=(1/N)Σi=1N(xix¯)(yyy¯)
{}corr({xi},{yi})