Generierung korrelierter binomialer Zufallsvariablen

Ich habe mich gefragt, ob es nach einem linearen Transformationsansatz möglich sein könnte, korrelierte binomische Zufallsvariablen zu generieren.

Unten habe ich etwas Einfaches in R ausprobiert und es ergibt eine gewisse Korrelation. Aber ich habe mich gefragt, ob es einen prinzipiellen Weg gibt, dies zu tun?

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)

Binomialvariablen werden normalerweise durch Summieren unabhängiger Bernoulli-Variablen erstellt. Mal sehen, ob wir mit einem Paar korrelierter Bernoulli-Variablen und dasselbe tun können.$(X,Y)$

Angenommen, ist eine Bernoulli -Variable ( und ) und ist eine Bernoulli -Variable. Um ihre gemeinsame Verteilung festzulegen, müssen wir alle vier Kombinationen von Ergebnissen angeben. Wenn wir schreiben wir den Rest leicht aus den Axiomen der Wahrscheinlichkeit herausfinden:( p ) Pr ( X = 1 ) = p Pr ( X = 0 ) = 1 - p Y ( q ) Pr ( ( X , Y ) = ( 0 , 0 ) ) = a , Pr ( ( X , Y ) = ( 1 , 0 ) ) = 1 - $X$$(p)$$\Pr(X=1)=p$$\Pr(X=0)=1-p$$Y$$(q)$

Das Einfügen in die Formel für den Korrelationskoeffizienten und das Lösen ergibt$\rho$

Vorausgesetzt, alle vier Wahrscheinlichkeiten sind nicht negativ, ergibt dies eine gültige gemeinsame Verteilung - und diese Lösung parametrisiert alle bivariaten Bernoulli-Verteilungen. (Wenn , gibt es eine Lösung für alle mathematisch bedeutsamen Korrelationen zwischen und ) Wenn wir dieser Variablen summieren , bleibt die Korrelation gleich - aber jetzt sind die Randverteilungen Binomial und Binomial , wie gewünscht.$p=q$$-1$$1$$n$$(n,p)$$(n,q)$

Beispiel

Es sei , , , und wir möchten, dass die Korrelation . Die Lösung ist (und die anderen Wahrscheinlichkeiten sind etwa , und ). Hier ist eine Darstellung von Realisierungen aus der gemeinsamen Verteilung:$n=10$$p=1/3$$q=3/4$$\rho=-4/5$$(1)$$a=0.00336735$$0.247$$0.663$$0.087$$1000$

Die roten Linien geben das Mittel der Stichprobe an und die gepunktete Linie ist die Regressionslinie. Sie sind alle nahe an ihren beabsichtigten Werten. Die Punkte in diesem Bild wurden zufällig gezittert, um die Überlappungen aufzulösen: Schließlich erzeugen Binomialverteilungen nur ganzzahlige Werte, sodass es zu einer starken Überzeichnung kommt.

Eine Möglichkeit, diese Variablen zu generieren, besteht darin, mal aus mit den gewählten Wahrscheinlichkeiten abzutasten und dann jede in , jede in , jede umzuwandeln in und jeweils in . Summiere die Ergebnisse (als Vektoren), um eine Realisierung von .{ 1 , 2 , 3 , 4 } 1 ( 0 , 0 ) 2 ( 1 , 0 ) $n$$\{1,2,3,4\}$$1$$(0,0)$$2$$(1,0)$$3$4 ( 1 , 1 ) ( X , Y $(0,1)$$4$$(1,1)$$(X,Y)$

Code

Hier ist eine RImplementierung.

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