Generuj pary liczb losowych równomiernie rozmieszczonych i skorelowanych

14

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?

Onturenio
źródło
6
Ściśle powiązane: stats.stackexchange.com/questions/30526 . Chcesz również sprawdzić tag kopuły - wystarczy kliknąć link tutaj. Szybka i brudna technika polega na tym, aby był jednolity [ 0 , 1 ], a Y = X, gdy X α, a Y = 1 + α - X w przeciwnym razie. Korelacja wynosi ρ = 2 ( α - 1 ) 3 + 1 , skąd α = 1 -X[0,1]Y=XXαY=1+αXρ=2(α1)3+1α=1((1ρ)/2)1/3Zrób sztuczkę. Ale kopuły dadzą ci większą kontrolę ...
whuber
Dziękuję za komentarz, ale tak, myślę, że ta metoda jest naprawdę „brudna”
Onturenio
1
Mam nadzieję, że widząc to podejście, rozpoznasz, że możesz (i powinieneś) podać dodatkowe kryteria dotyczące właściwości twoich par liczb losowych. Jeśli jest to „brudne”, to co dokładnie jest nie tak z rozwiązaniem? Powiedz nam, abyśmy mogli udzielić bardziej odpowiednich odpowiedzi w Twojej sytuacji.
whuber
Odpowiedzi na to pytanie udzielono przypadkowo w odpowiedzi na ściśle powiązane pytanie: jak wygenerować pary RV z regresją liniową. Ponieważ nachylenie regresji liniowej jest powiązane w łatwo obliczalny sposób ze współczynnikiem korelacji i można wygenerować wszystkie możliwe nachylenia, daje to sposób na wytworzenie dokładnie tego, co chcesz. Zobacz stats.stackexchange.com/questions/257779/… .
whuber
1
Zobacz także stats.stackexchange.com/questions/31771 , który odpowiada uogólnieniu na trzy losowe mundury.
whuber

Odpowiedzi:

16

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 ) =U1U2F1F2Fi(Ui)=Uija=1,2) Tak więc rho Spearmana i współczynnik korelacji Pearsona są równe (wersje przykładowe mogą się jednak różnić).

ρS(U1,U2)=corr(F1(U1),F2(U2))=corr(U1,U2).

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 , XX1,X2ρ Ułatwia to generowanie losowych zmiennych, które mają pożądaną wartość rho Spearmana.

ρS(X1,X2)=6πarcsin(ρ2).

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.ρ

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

  1. Oblicz .ρ=2sin(rπ/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.6n=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 iU1U2U1 . U2wprowadź opis zdjęcia tutaj

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ównyr .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.6n

## 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")

wprowadź opis zdjęcia tutaj wprowadź opis zdjęcia tutaj wprowadź opis zdjęcia tutaj

QuantIbex
źródło
3
Ogólny sposób generowania skorelowanych rozkładów wielowymiarowych z danymi rozkładami krańcowymi nazywa się kopulą .
whuber
@ whuber użycie kopuły pozwala określić strukturę zależności między zmiennymi losowymi. Problem polega na tym, że na korelację (Osoba) wpływają zarówno struktura zależności, jak i marginesy. Zatem każdy wybór marginesów będzie wymagał odpowiedniego wyboru parametrów kopuły, nie wspominając już o tym, że niektórych poziomów korelacji po prostu nie można osiągnąć dla danych marginesów (np. Patrz tutaj ). Jeśli znasz metodę, która pozwala „kontrolować” poziom korelacji dla dowolnego wyboru marginesów, chciałbym o tym wiedzieć.
QuantIbex,
Dzięki @QuantIbex. Ale nie rozumiem, dlaczego „a) sugeruje, że współczynnik korelacji rho i (Pearsona) Spearmana dla zmiennych losowych ze standardowymi jednolitymi marginesami jest w przybliżeniu równy w dużej próbce”
Onturenio
2
[1,1]
1
@Quantibex Pozwoliłem sobie dodać zdanie, które wskazuje, że twoja gen.gauss.copfunkcja będzie działać dla więcej niż dwóch zmiennych z (trywialną) poprawką. Jeśli nie podoba ci się ten dodatek lub chcesz go inaczej ująć, cofnij lub zmień w razie potrzeby.
Glen_b
0

u1U(0,1)u1w1U(0,1)I=1u1w2U(0,1)I=0u1U(0,1)u2

E(u1u2)=E[Iw1+(1I)w2][Iw1+(1I)w3]

I(I1)=0I2=I(1I)2=(1I)I01Iw

E(u1u2)=E(I)E(w12)+E(1I)E(w2)E(w3) =pE(w12)+(1p)/4

V(w1)=1/12E(w12)=1/3E(u1u2)=p/12+1/4cov(u1u2)=p/12. Since V(u1)=V(u2)=1/12, we get finally that cor(u1,u2)=p.

Neal Oden
źródło
0

Here is one easy method for positive correlation: Let (u1,u2)=Iw1+(1I)(w2,w3), where w1,w2, and w3 are independent U(0,1) and I is Bernoulli(p). u1 and u2 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 (u1,u2)=I(w1,1w1)+(1I)(w2,w3), and the correlation will be p.

Neal Oden
źródło
Can you add a short proof of why this works?
The Laconic
if your want to be computationally efficient, u1=w1 also produces the same correlation (both positive and negative cases)
Anvit