Erzeugen Sie Paare von Zufallszahlen, die gleichmäßig verteilt und korreliert sind

Ich möchte Zufallszahlenpaare mit einer bestimmten Korrelation erzeugen. Der übliche Ansatz, eine Linearkombination zweier Normalvariablen zu verwenden, ist hier jedoch nicht gültig, da eine Linearkombination gleichförmiger Variablen keine gleichmäßig verteilte Variable mehr ist. Ich brauche die beiden Variablen, um einheitlich zu sein.

Irgendeine Idee, wie Paare von einheitlichen Variablen mit einer gegebenen Korrelation erzeugt werden können?

Mir ist keine universelle Methode bekannt, um korrelierte Zufallsvariablen mit einer gegebenen Randverteilung zu erzeugen. Ich werde daher eine Ad-hoc-Methode vorschlagen, um Paare gleichmäßig verteilter Zufallsvariablen mit einer gegebenen (Pearson-) Korrelation zu generieren. Ohne Einschränkung der Allgemeinheit gehe ich davon aus, dass die gewünschte Randverteilung einheitlich ist (dh der Träger ist ).$[0, 1]$

Der vorgeschlagene Ansatz beruht auf Folgendem:
a) Für standardmäßige einheitliche Zufallsvariablen und mit den jeweiligen Verteilungsfunktionen und gilt , für . Somit kann durch Definition Spearman-rho ist Daher sind Spearmans Rho und Pearsons Korrelationskoeffizient gleich (Beispielversionen können sich jedoch unterscheiden).$U_1$$U_2$$F_1$$F_2$$F_i(U_i) = U_i$$i = 1, 2$

b) Wenn Zufallsvariablen mit stetigen Rändern und Gaußsche Kopula mit (Pearson) -Korrelationskoeffizient , dann ist Spearmans rho Dies macht es einfach, Zufallsvariablen zu generieren, die einen gewünschten Wert von Spearmans Rho haben.$X_1, X_2$$\rho$

Der Ansatz besteht darin, Daten aus der Gaußschen Kopula mit einem geeigneten Korrelationskoeffizienten zu generieren, so dass das Spearman-Rho der gewünschten Korrelation für die einheitlichen Zufallsvariablen entspricht.$\rho$

Simulationsalgorithmus
Es sei der gewünschte Korrelationsgrad und die Anzahl der zu erzeugenden Paare. Der Algorithmus ist:$r$$n$

1. Berechne .$\rho = 2\sin (r \pi/6)$
2. Generieren Sie ein Paar zufälliger Variablen aus der Gaußschen Copula (z. B. mit diesem Ansatz )
3. Wiederholen Sie Schritt 2 Mal.$n$

Beispiel
Der folgende Code ist ein Beispiel für die Implementierung dieses Algorithmus unter Verwendung von R mit einer Zielkorrelation von und Paaren.$r = 0.6$$n = 500$

## Initialization and parameters 
set.seed(123)
r <- 0.6                            # Target (Spearman) correlation
n <- 500                            # Number of samples

## Functions
gen.gauss.cop <- function(r, n){
    rho <- 2 * sin(r * pi/6)        # Pearson correlation
    P <- toeplitz(c(1, rho))        # Correlation matrix
    d <- nrow(P)                    # Dimension
    ## Generate sample
    U <- pnorm(matrix(rnorm(n*d), ncol = d) %*% chol(P))
    return(U)
}

## Data generation and visualization
U <- gen.gauss.cop(r = r, n = n)
pairs(U, diag.panel = function(x){
          h <- hist(x, plot = FALSE)
          rect(head(h$breaks, -1), 0, tail(h$breaks, -1), h$counts/max(h$counts))})

In der folgenden Abbildung zeigen die Diagonaldiagramme Histogramme der Variablen und und die nicht diagonalen Diagramme Streudiagramme von und . $U_1$$U_2$$U_1$$U_2$

Konstruktiv haben die Zufallsvariablen einheitliche Ränder und einen Korrelationskoeffizienten (nahe) . Aufgrund der Auswirkung der Abtastung ist der Korrelationskoeffizient der simulierten Daten jedoch nicht genau gleich .$r$$r$

cor(U)[1, 2]
# [1] 0.5337697

Beachten Sie, dass die gen.gauss.copFunktion mit mehr als zwei Variablen arbeiten sollte, indem Sie einfach eine größere Korrelationsmatrix angeben.

Simulationsstudie
Die folgende Simulationsstudie, die für die Zielkorrelation wiederholt legt nahe, dass die Verteilung des Korrelationskoeffizienten mit zunehmender Stichprobengröße gegen die gewünschte Korrelation konvergiert .$r= -0.5, 0.1, 0.6$$n$

## Simulation
set.seed(921)
r <- 0.6                                                # Target correlation
n <- c(10, 50, 100, 500, 1000, 5000); names(n) <- n     # Number of samples
S <- 1000                                               # Number of simulations

res <- sapply(n,
              function(n, r, S){
                   replicate(S, cor(gen.gauss.cop(r, n))[1, 2])
               }, 
               r = r, S = S)
boxplot(res, xlab = "Sample size", ylab = "Correlation")
abline(h = r, col = "red")

Intuitiv ist U ( 0 , 1 ) , da u 1 equals w 1 [welche U ( 0 , 1 ) ] , wenn I = 1 , und u 1 equals w 2 [welche U ( 0 , 1 ) ] , wenn I = 0 , so u 1 ist U ( 0 , 1 $u_1$$U(0,1)$$u_1$$w_1$$U(0,1)$$I = 1$$u_1$$w_2$$U(0,1)$$I = 0$$u_1$$U(0,1)$in beiden Fällen. Das selbe für . Zur Korrelation:$u_2$

$E(u_1 u_2) = E[I w_1 + (1-I) w_2][I w_1 + (1-I) w_3]$

$I(I-1)=0$$I^2=I$$(1-I)^2=(1-I)$$I$$0$$1$$I$$w$

$E(u_1 u_2) = E(I)E(w_1^2) + E(1-I)E(w_2)E(w_3)$ $=pE(w_1^2)+(1-p)/4$

$V(w_1)=1/12$$E(w_1^2)=1/3$$E(u_1 u_2) = p/12 + 1/4$$cov(u_1 u_2) = p/12$$V(u_1)=V(u_2)=1/12$$cor(u_1, u_2) = p$

$(u_1, u_2) = Iw_1 + (1-I) (w_2, w_3)$$w_1, w_2,$$w_3$$U(0,1)$$I$$p$$u_1$$u_2$$U(0,1)$$p$$k$

$(u_1, u_2) = I(w_1, 1-w_1) + (1-I)(w_2, w_3)$, and the correlation will be $-p$.