Generuj pary liczb losowych równomiernie rozmieszczonych i skorelowanych

Chciałbym wygenerować pary liczb losowych z pewną korelacją. Jednak zwykłe podejście polegające na stosowaniu kombinacji liniowej dwóch zmiennych normalnych nie jest tutaj poprawne, ponieważ kombinacja liniowa zmiennych jednolitych nie jest już zmienną równomiernie rozłożoną. Potrzebuję dwóch zmiennych, aby były jednolite.

Masz pomysł, jak wygenerować pary zmiennych jednorodnych o zadanej korelacji?

Nie znam uniwersalnej metody generowania skorelowanych zmiennych losowych z dowolnymi rozkładami krańcowymi. Tak więc zaproponuję metodę ad hoc do generowania par równomiernie rozmieszczonych zmiennych losowych o zadanej korelacji (Pearsona). Bez utraty ogólności, zakładam, że pożądany rozkład krańcowy jest standardowy jednolity (tj. Wsparcie wynosi ).$[0, 1]$

Proponowane podejście opiera się na:
a) Dla standardowych jednorodnych zmiennych losowych i U 2 z odpowiednimi funkcjami rozkładu F 1 i F 2 , mamy F i ( U i ) = U i , dla i = 1 , 2 . Zatem z definicji rho Spearmana wynosi ρ S ( U 1 , U 2 ) $U_1$$U_2$$F_1$$F_2$$F_i(U_i) = U_i$$i = 1, 2$ Tak więc rho Spearmana i współczynnik korelacji Pearsona są równe (wersje przykładowe mogą się jednak różnić).

b) Jeśli są zmiennymi losowymi z ciągłymi marginesami i kopulą Gaussa ze współczynnikiem korelacji (Pearsona) ρ , wówczas rho Spearmana wynosi ρ S ( X 1 , $X_1, X_2$$\rho$ Ułatwia to generowanie losowych zmiennych, które mają pożądaną wartość rho Spearmana.

Podejście polega na generowaniu danych z kopuły Gaussa z odpowiednim współczynnikiem korelacji tak aby rho Spearmana odpowiadało pożądanej korelacji dla jednolitych zmiennych losowych.$\rho$

Algorytm symulacji
Niech $r$ oznacza pożądany poziom korelacji, a liczbę generowanych par. Algorytm to:$n$

1. Oblicz .$\rho = 2\sin (r \pi/6)$
2. Wygeneruj parę zmiennych losowych z kopuły Gaussa (np. Przy takim podejściu )
3. Powtórz krok 2 razy.$n$

Przykład
Poniższy kod jest przykładem implementacji tego algorytmu przy użyciu R z korelacją docelową i n = 500 par.$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))})

Na poniższym rysunku wykresy ukośne pokazują histogramy zmiennych i U 2 , a wykresy nie przekątne pokazują wykresy rozproszenia U 1 i$U_1$$U_2$$U_1$ . $U_2$

Konstruując, zmienne losowe mają jednolite marginesy i współczynnik korelacji (bliski) . Ale ze względu na efekt próbkowania współczynnik korelacji symulowanych danych nie jest dokładnie równy$r$ .$r$

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

Zauważ, że gen.gauss.copfunkcja powinna działać z więcej niż dwiema zmiennymi, po prostu określając większą macierz korelacji.

Badanie symulacyjne
Poniższe badanie symulacyjne powtórzone dla korelacji docelowej sugeruje, że rozkład współczynnika korelacji zbliża się do pożądanej korelacji wraz ze wzrostem wielkości próby n .$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")

$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)$$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$. Since $V(u_1)=V(u_2)=1/12$, we get finally that $cor(u_1, u_2) = p$.

Here is one easy method for positive correlation: Let $(u_1, u_2) = Iw_1 + (1-I) (w_2, w_3)$, where $w_1, w_2,$ and $w_3$ are independent $U(0,1)$ and $I$ is Bernoulli($p$). $u_1$ and $u_2$ will then have $U(0,1)$ distributions with correlation $p$. This extends immediately to $k$-tuples of uniforms with compound symmetric variance matrix.

If you want pairs with negative correlation, use $(u_1, u_2) = I(w_1, 1-w_1) + (1-I)(w_2, w_3)$, and the correlation will be $-p$.