Verwenden der Gleichverteilung zur Erzeugung korrelierter Zufallsstichproben in R.

[Bei den letzten Fragen habe ich mich mit der Erzeugung von Zufallsvektoren in R befasst und wollte diese "Forschung" als unabhängige Frage und Antwort zu einem bestimmten Punkt teilen.]

Erzeugen von Zufallsdaten mit Korrelation kann unter Verwendung der Cholesky - Zerlegung der Korrelationsmatrix durchgeführt wird hier , wie auf dem Stand der Beiträge reflektierte hier und hier . $C = LL^{T}$

Die Frage, die ich ansprechen möchte, ist, wie die Gleichverteilung verwendet werden kann, um korrelierte Zufallszahlen aus verschiedenen Randverteilungen in R zu generieren .

r correlation sampling random-variable random-generation Antoni Parellada
quelle

Sie scheinen die Gaußsche Kopula wiederentdeckt zu haben, z. B. siehe verwandte Frage hier . Es gibt viele andere Copulas, die häufig verwendet werden, aber der Gaußsche ist sehr praktisch und kann für einige Situationen durchaus geeignet sein.

Glen_b -State Monica

Antworten:

Da ist die Frage

"Wie verwende ich die Gleichverteilung, um korrelierte Zufallszahlen aus verschiedenen Randverteilungen in zu generieren? " $\mathbb{R}$

und nicht nur normale Zufallsvariablen, die obige Antwort erzeugt keine Simulationen mit der beabsichtigten Korrelation für ein beliebiges Paar von Randverteilungen in . $\mathbb{R}$

Der Grund ist , dass für die meisten cdfs und , wenn wobei das normale Standard-cdf bezeichnet. $G_X$ $G_Y$

cor (X, Y) \neq cor (G_{X}^{- 1} (Φ (X), G_{Y}^{- 1} (Φ (Y)),

$\text{cor}(X,Y)\ne\text{cor}(G_X^{-1}(\Phi(X),G_Y^{-1}(\Phi(Y)),$

(X, Y) \sim N_{2} (0, Σ),

$(X,Y)\sim\text{N}_2(0,\Sigma),$

Φ

$\Phi$

Nämlich ist hier ein Gegenbeispiel mit einem Exp (1) und eine Gamma (.2,1) als mein Paar Randverteilungen in . $\mathbb{R}$

library(mvtnorm)
#correlated normals with correlation 0.7
x=rmvnorm(1e4,mean=c(0,0),sigma=matrix(c(1,.7,.7,1),ncol=2),meth="chol")
cor(x[,1],x[,2])
  [1] 0.704503
y=pnorm(x) #correlated uniforms
cor(y[,1],y[,2])
  [1] 0.6860069
#correlated Exp(1) and Ga(.2,1)
cor(-log(1-y[,1]),qgamma(y[,2],shape=.2))
  [1] 0.5840085

Ein weiteres offensichtliches Gegenbeispiel ist, wenn das Cauchy-cdf ist. In diesem Fall ist die Korrelation nicht definiert. $G_X$

Um ein breiteres Bild zu geben, hier ein R-Code, bei dem sowohl als auch beliebig sind: $G_X$ $G_Y$

etacor=function(rho=0,nsim=1e4,fx=qnorm,fy=qnorm){
  #generate a bivariate correlated normal sample
  x1=rnorm(nsim);x2=rnorm(nsim)
  if (length(rho)==1){
    y=pnorm(cbind(x1,rho*x1+sqrt((1-rho^2))*x2))
    return(cor(fx(y[,1]),fy(y[,2])))
    }
  coeur=rho
  rho2=sqrt(1-rho^2)
  for (t in 1:length(rho)){
     y=pnorm(cbind(x1,rho[t]*x1+rho2[t]*x2))
     coeur[t]=cor(fx(y[,1]),fy(y[,2]))}
  return(coeur)
  }

Geben Sie hier die Bildbeschreibung ein

Das Herumspielen mit verschiedenen cdfs führte mich dazu, diesen Sonderfall einer -Verteilung für und einer logarithmischen Normalverteilung für herauszustellen : $\chi^2_3$ $G_X$ $G_Y$

rhos=seq(-1,1,by=.01)
trancor=etacor(rho=rhos,fx=function(x){qchisq(x,df=3)},fy=qlnorm)
plot(rhos,trancor,ty="l",ylim=c(-1,1))
abline(a=0,b=1,lty=2)

Dies zeigt, wie weit die Korrelation von der Diagonale entfernt sein kann.

$G_X$ $G_Y$ $\text{cor}(X,Y)$ $(-1,1)$

Xi'an
quelle

Fantastisch! Ty! Gibt es eine Möglichkeit, ein ungefähres Segment zu finden, in dem die Abweichung nicht markiert ist, wie es bei Normalen der Fall zu sein scheint, um für praktische Anwendungen noch angemessen zu sein?

Antoni Parellada

Ich habe das correlatePaket geschrieben. Die Leute sagten, es sei vielversprechend (verdient eine Veröffentlichung im Journal of Statistical Software), aber ich habe das Papier nie dafür geschrieben, weil ich mich entschieden habe, keine akademische Karriere zu verfolgen.

Ich glaube, das nicht gewartete correlatePaket ist immer noch auf CRAN.

Wenn Sie es installieren, können Sie Folgendes tun:

require('correlate')
a <- rnorm(100)
b <- runif(100)
newdata <- correlate(cbind(a,b),0.5)

Das Ergebnis ist, dass neue Daten eine Korrelation von 0,5 haben, ohne die univariaten Verteilungen von aund zu ändern b(dieselben Werte sind vorhanden, sie werden nur verschoben, bis die multivariate 0,5-Korrelation erreicht ist.

Ich werde hier auf Fragen antworten, entschuldige die fehlende Dokumentation.

PascalVKooten
quelle

Bravo, das ist die perfekte Antwort! Haben Sie eine Möglichkeit, Werte der Korrelation zu erkennen, die nicht zu erreichen sind?

Xi'an

@ Xi'an Es gibt einige Unmöglichkeiten, wie wenige Datenpunkte und eine wirklich spezifische Korrelation, die einfach nicht erreicht werden kann. zB nur mit 3 gepaarten Werten.

PascalVKooten

Beachten Sie auch, dass es für mehr als 2 Variablen möglich ist, z. B. für 3 Variablen können Sie eine 3x3-Korrelationsmatrix definieren, 4 Variablen eine 4x4.

PascalVKooten

Im Allgemeinen funktioniert es, solange Sie nicht das Unmögliche wollen, aber bevor Sie ernsthaft damit arbeiten, wird empfohlen, einige Testläufe durchzuführen.

PascalVKooten

Leute, die daran interessiert waren, verwendeten Einkommensdaten; Ladungen von Nullen und eine Gaußsche Verteilung für Einkommen ungleich Null.

PascalVKooten

Generieren Sie zwei Stichproben korrelierter Daten aus einer normalen Standard-Zufallsverteilung nach einer vorgegebenen Korrelation .

Als Beispiel wählen wir eine Korrelation r = 0,7 und codieren eine Korrelationsmatrix wie:

(C <- matrix(c(1,0.7,0.7,1), nrow = 2)) [,1] [,2] [1,] 1.0 0.7 [2,] 0.7 1.0

Wir können mvtnormjetzt diese beiden Stichproben als bivariaten Zufallsvektor generieren:

set.seed(0)

SN <- rmvnorm(mean = c(0,0), sig = C, n = 1e5) $N(0, 1)$ cor(SN[,1],SN[,2])= 0.6996197 ~ 0.7

X1 <- SN[,1]; X2 <- SN[,2]

Hier ist die Darstellung mit der überlappenden Regressionslinie:
$U(0, 1)$

U <- pnorm(SN)pnormSN $erf(SN)$ $\Phi(SN)$ cor(U[,1], U[,2]) = 0.6816123 ~ 0.7

Wieder können wir den Vektor zerlegen U1 <- U[,1]; U2 <- U[,2]und ein Streudiagramm mit Randverteilungen an den Rändern erzeugen, das ihre einheitliche Natur deutlich zeigt:
Anwenden der inversen Transformationsverfahren Abtasten hier , um schließlich den Bivektor von gleich korrelierten Punkte zu erhalten , die zu je nachdem , welche Verteilung Familie wir uns vorgenommen haben zu reproduzieren.

Von hier aus können wir nur zwei Vektoren erzeugen, die normal und mit gleichen oder unterschiedlichen Varianzen verteilt sind . Zum Beispiel: Y1 <- qnorm(U1, mean = 8,sd = 10)und Y2 <- qnorm(U2, mean = -5, sd = 4), wodurch die gewünschte Korrelation aufrechterhalten wird , cor(Y1,Y2) = 0.6996197 ~ 0.7.

U1 $t$ U2 $\lambda$ Z1 <- qt(U1, df = 3)Z2 <- qexp(U2, rate = 1)cor(Z1,Z2) [1] 0.5941299 < 0.7

Hier ist ein Beispiel für Code für den gesamten Prozess und normale Ränder:

Cor_samples <- function(r, n, mean1, mean2, sd1, sd2){
C <- matrix(c(1,r,r,1), nrow = 2)
require(mvtnorm)
SN <- rmvnorm(mean = c(0,0), sig = C, n = n)
U <- pnorm(SN)
U1 <- U[,1]
U2 <- U[,2]

 Y1 <<- qnorm(U1, mean = mean1,sd = sd1) 
 Y2 <<- qnorm(U2, mean = mean2,sd = sd2) 

sample_measures <<- as.data.frame(c(mean(Y1), mean(Y2), sd(Y1), sd(Y2), cor(Y1,Y2)), names<-c("mean Y1", "mean Y2", "SD Y1", "SD Y2", "Cor(Y1,Y2)"))
sample_measures
}

Zum Vergleich habe ich eine Funktion zusammengestellt, die auf der Cholesky-Zerlegung basiert:

Cholesky_samples <- function(r, n, mean1, mean2, sd1, sd2){
C <- matrix(c(1,r,r,1), nrow = 2)
L <- chol(C)
X1 <- rnorm(n)
X2 <- rnorm(n)
X <- rbind(X1,X2)

Y <- t(L)%*%X
Y1 <- Y[1,]
Y2 <- Y[2,]

N_1 <<- Y[1,] * sd1 + mean1
N_2 <<- Y[2,] * sd2 + mean2

sample_measures <<- as.data.frame(c(mean(N_1), mean(N_2), sd(N_1), sd(N_2), cor(N_1, N_2)), 
                  names<-c("mean N_1", "mean N_2", "SD N_1", "SD N_2","cor(N_1,N_2)"))
sample_measures
}

$r=0.7$ $N(97,23)$ $N(32,8)$ set.seed(99)

Verwenden der Uniform:

cor_samples(0.7, 1000, 97, 32, 23, 8)
           c(mean(Y1), mean(Y2), sd(Y1), sd(Y2), cor(Y1, Y2))
mean Y1                                            96.5298821
mean Y2                                            32.1548306
SD Y1                                              22.8669448
SD Y2                                               8.1150780
cor(Y1,Y2)                                          0.7061308

und Verwenden des Cholesky:

Cholesky_samples(0.7, 1000, 97, 32, 23, 8)
             c(mean(N_1), mean(N_2), sd(N_1), sd(N_2), cor(N_1, N_2))
mean N_1                                                   96.4457504
mean N_2                                                   31.9979675
SD N_1                                                     23.5255419
SD N_2                                                      8.1459100
cor(N_1,N_2)                                                0.7282176

Antoni Parellada
quelle

F^{- 1} (X)

$F^{-1}(X)$

f (F^{- 1} (X))

$f(F^{-1}(X))$

Ich habe die Funktion am Ende der Antwort geändert, um die Korrelation der berechneten Stichproben einzuschließen, um sie mit der eingesteckten Zahl zu vergleichen, und sie scheinen übereinzustimmen.

Antoni Parellada

Ob es Probleme mit der praktischen Anwendung gibt, hängt von der praktischen Anwendung ab. Für einige Dinge ist das in Ordnung. Da die Transformationen monoton sind, werden nichtparametrische Korrelationen wie Spearmans Rho und Kendalls Tau nicht geändert.

Glen_b -Reinstate Monica