Generowanie skorelowanych zmiennych losowych dwumianowych

Zastanawiałem się, czy możliwe jest wygenerowanie skorelowanych losowych zmiennych dwumianowych po zastosowaniu transformacji liniowej?

Poniżej wypróbowałem coś prostego w R i daje to pewną korelację. Ale zastanawiałem się, czy istnieje jakiś zasadny sposób, aby to zrobić?

X1 = rbinom(1e4, 6, .5) ; X2 = rbinom(1e4, 6, .5) ;  X3 = rbinom(1e4, 6, .5) ; a = .5

Y1 = X1 + (a*X2) ; Y2 = X2 + (a*X3) ## Y1 and Y2 are supposed to be correlated

cor(Y1, Y2)

— rnorouzian
źródło

Y_{1}

$Y_1$ i mogą być skorelowane, ale nie będą już dwumianowe. Przykład: a następnie stąd nie może być losowymi zmiennymi dwumianowymi. Sugerowałbym przyjrzenie się rozkładowi wielomianowemu.

Y_{2}

$Y_2$

X_{1} = X_{2} = 1

$X_1 = X_2 = 1$

Y_{1} = 1.5

$Y_1 = 1.5$

Y_{i}

$Y_i$

— knrumsey - Przywróć Monikę

Krótka odpowiedź na pytanie polega na szukaniu słowa kluczowego copula, które pomaga w generowaniu zmiennych zależnych o ustalonych marginesach.

— Xi'an

Zmienne dwumianowe są zwykle tworzone przez zsumowanie niezależnych zmiennych Bernoulliego. Zobaczmy, czy możemy zacząć od pary skorelowanych zmiennych Bernoulliego i zrobić to samo. $(X,Y)$

Załóżmy, że jest zmienną Bernoulliego (to znaczy a ), a jest zmienną Bernoulliego . Aby określić ich wspólną dystrybucję, musimy określić wszystkie cztery kombinacje wyników. Pisząc możemy łatwo ustalić resztę z aksjomatów prawdopodobieństwa: $X$ $(p)$ $\Pr(X=1)=p$ $\Pr(X=0)=1-p$ $Y$ $(q)$

Par ((X, Y) = (0, 0)) = za,

$\Pr((X,Y)=(0,0))=a,$

Par ((X, Y) = (1, 0)) = 1 - q - za, Par ((X, Y) = (0, 1)) = 1 - p - za, Par ((X, Y) = (1, 1)) = za + p + q - 1.

$\Pr((X,Y)=(1,0))=1-q-a, \\\Pr((X,Y)=(0,1))=1-p-a, \\\Pr((X,Y)=(1,1))=a+p+q-1.$

Podłączenie tego do wzoru dla współczynnika korelacji i rozwiązanie daje $\rho$

\begin{matrix} (1) & za = (1 - p) (1 - q) + ρ \sqrt{p q (1 - p) (1 - q)} . \end{matrix}

$a = (1-p)(1-q) + \rho\sqrt{{pq}{(1-p)(1-q)}}.\tag{1}$

Pod warunkiem, że wszystkie cztery prawdopodobieństwa nie są ujemne, da to prawidłowy wspólny rozkład - a to rozwiązanie parametryzuje wszystkie dwuwymiarowe rozkłady Bernoulliego. (Gdy , istnieje rozwiązanie dla wszystkich matematycznie znaczących korelacji między a ) Gdy sumujemy tych zmiennych, korelacja pozostaje taka sama - ale teraz rozkłady brzeżne są dwumianowe i Dwumianowy , zgodnie z potrzebami. $p=q$ $-1$ $1$ $n$ $(n,p)$ $(n,q)$

Przykład

Niech , , , i chcielibyśmy, aby korelacja wynosiła . Rozwiązaniem dla jest (a inne prawdopodobieństwa wynoszą około , i ). Oto wykres realizacji ze wspólnej dystrybucji: $n=10$ $p=1/3$ $q=3/4$ $\rho=-4/5$ $(1)$ $a=0.00336735$ $0.247$ $0.663$ $0.087$ $1000$

Czerwone linie wskazują średnie próbki, a linia przerywana to linia regresji. Wszystkie są zbliżone do zamierzonych wartości. Punkty zostały przypadkowo roztrzęsione na tym obrazie, aby rozwiązać nakładanie się: w końcu rozkłady dwumianowe produkują tylko wartości całkowite, więc będzie bardzo dużo wykreślania.

Jednym ze sposobów generowania tych zmiennych jest próbkowanie razy z z wybranymi prawdopodobieństwami, a następnie przeliczanie każdego na , każdego na , każde na i każdy na . Zsumuj wyniki (jako wektory), aby uzyskać jedną realizację . $n$ $\{1,2,3,4\}$ $1$ $(0,0)$ $2$ $(1,0)$ $3$ $(0,1)$ $4$ $(1,1)$ $(X,Y)$

Kod

Oto Rimplementacja.

#
# Compute Pr(0,0) from rho, p=Pr(X=1), and q=Pr(Y=1).
#
a <- function(rho, p, q) {
  rho * sqrt(p*q*(1-p)*(1-q)) + (1-p)*(1-q)
}
#
# Specify the parameters.
#
n <- 10
p <- 1/3
q <- 3/4
rho <- -4/5
#
# Compute the four probabilities for the joint distribution.
#
a.0 <- a(rho, p, q)
prob <- c(`(0,0)`=a.0, `(1,0)`=1-q-a.0, `(0,1)`=1-p-a.0, `(1,1)`=a.0+p+q-1)
if (min(prob) < 0) {
  print(prob)
  stop("Error: a probability is negative.")
}
#
# Illustrate generation of correlated Binomial variables.
#
set.seed(17)
n.sim <- 1000
u <- sample.int(4, n.sim * n, replace=TRUE, prob=prob)
y <- floor((u-1)/2)
x <- 1 - u %% 2
x <- colSums(matrix(x, nrow=n)) # Sum in groups of `n`
y <- colSums(matrix(y, nrow=n)) # Sum in groups of `n`
#
# Plot the empirical bivariate distribution.
#
plot(x+rnorm(length(x), sd=1/8), y+rnorm(length(y), sd=1/8),
     pch=19, cex=1/2, col="#00000010",
     xlab="X", ylab="Y",
     main=paste("Correlation is", signif(cor(x,y), 3)))
abline(v=mean(x), h=mean(y), col="Red")
abline(lm(y ~ x), lwd=2, lty=3)

— Whuber
źródło

Czy to podejście można rozszerzyć o generowanie dowolnej liczby zmiennych binarnych? - pasować do danej macierzy korelacji (lub maksymalnie ją dopasować)?

— ttnphns

@ttnphns Tak, ale opcje eksplodują: tabela prawdopodobieństwa musi być określona przez parametrów krańcowych, ograniczenie sumy do jedności i (dlatego) dodatkowe parametry . Oczywiście masz dużą swobodę wyboru (lub ograniczenia) tych parametrów, zgodnie z właściwościami wielu zmiennych, które chcesz utworzyć. Podobne podejście można również zastosować do wygenerowania skorelowanych zmiennych dwumianowych o różnych wartościach ich parametrów „ ”. Parvin: Uważam, że „kiedy sumujemy tych zmiennych” jednoznacznie wyjaśnia, co reprezentuje .

2^{k}

$2^k$

k

$k$

2^{k} - k - 1

$2^k-k-1$

n

$n$

n

$n$

n

$n$

— whuber

To niezły wynik. Wystarczy wybrać trochę w pierwszym zdaniu. Czy aby uzyskać dwumian z niezależnych zmiennych Bernoulliego, nie muszą mieć tego samego p? Nie ma to wpływu na to, co zrobiłeś, ponieważ jest to tylko motywacja do podejścia.

— Michael R. Chernick

@Michael Dziękuję - masz całkowitą rację. Innym powodem, dla którego nie ma ona wpływu na opisaną tutaj metodę, jest to, że metoda ta nadal wymaga sumowania zmiennych Bernoulliego za pomocą wspólnego parametru: parametr jest dla wszystkich zmiennych i dla wszystkich zmiennychAby post był dość krótki, nie przedstawiłem histogramów rozkładów krańcowych, aby wykazać, że rzeczywiście są one dwumianowe (ale tak naprawdę zrobiłem to w mojej oryginalnej analizie, aby upewnić się, że działają!).

p

$p$

X

$X$

q

$q$

Y

$Y$

— whuber

@whuber Ładne podejście! Czy możesz dać mi znać, jeśli jest jakiś papier, do którego mogę się odwoływać?

— T Nick