Génération de variables aléatoires binomiales corrélées

Je me demandais s'il serait possible de générer des variables binomiales aléatoires corrélées en suivant une approche de transformation linéaire?

Ci-dessous, j'ai essayé quelque chose de simple en R et cela produit une certaine corrélation. Mais je me demandais s'il y avait un moyen de principe de le faire?

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)

Les variables binomiales sont généralement créées en additionnant des variables de Bernoulli indépendantes. Voyons si nous pouvons commencer avec une paire de variables de Bernoulli corrélées et faire la même chose.$(X,Y)$

Supposons que est une variable de Bernoulli ( p ) (c'est-à-dire Pr ( X = 1 ) = p et Pr ( X = 0 ) = 1 - p ) et Y est une variable de Bernoulli ( q ) . Pour déterminer leur distribution conjointe, nous devons spécifier les quatre combinaisons de résultats. Écriture Pr ( ( X , Y ) = ( 0 , 0 ) ) $X$$(p)$$\Pr(X=1)=p$$\Pr(X=0)=1-p$$Y$$(q)$ a

Le brancher dans la formule du coefficient de corrélation et résoudre donne a = ( 1 - p ) ( 1 $\rho$

À condition que les quatre probabilités ne soient pas négatives, cela donnera une distribution conjointe valide - et cette solution paramètre toutes les distributions bivariées de Bernoulli. (Lorsque , il existe une solution pour toutes les corrélations mathématiquement significatives entre - 1 et 1. ) Lorsque nous additionnons n de ces variables, la corrélation reste la même - mais maintenant les distributions marginales sont binomiales ( n , p ) et Binôme ( n , q $p=q$$-1$$1$$n$$(n,p)$$(n,q)$ , comme souhaité.

Exemple

Soit , p = 1 / 3 , q = 3 / 4 , et nous voudrions que la corrélation d'être ρ = - 4 / 5 . La solution de ( 1 ) est a = 0,00336735 (et les autres probabilités se situent autour de 0,247 , 0,663 et $n=10$$p=1/3$$q=3/4$$\rho=-4/5$$(1)$$a=0.00336735$$0.247$$0.663$$0.087$ ). Voici un tracé de réalisations de la distribution conjointe:$1000$

Les lignes rouges indiquent les moyennes de l'échantillon et la ligne pointillée est la ligne de régression. Ils sont tous proches de leurs valeurs prévues. Les points ont été agités au hasard dans cette image pour résoudre les chevauchements: après tout, les distributions binomiales ne produisent que des valeurs intégrales, il y aura donc une grande quantité de surplotage.

$n$$\{1,2,3,4\}$$1$$(0,0)$$2$$(1,0)$$3$$(0,1)$$4$$(1,1)$$(X,Y)$

Code

Voici une Rimplémentation.

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