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

21

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)
rnorouzian
la source
2
Y1 etY2 peuvent être corrélés, mais ils ne seront plus binomiaux. Exemple,X1=X2=1 puisY1=1.5 où leYi ne peut pas être des variables aléatoires binomiales. Je vous suggère de vous pencher sur la distribution multinomiale.
knrumsey
1
La réponse courte à la question est de rechercher le mot-clé copula, qui aide à générer des variables dépendantes avec des marges fixes.
Xi'an

Réponses:

32

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)=pPr(X=0)=1pY(q) a ,

Pr((X,Y)=(0,0))=a,
nous pouvons facilement comprendre le reste à partir des axiomes de probabilité:
Pr((X,Y)=(1,0))=1qa,Pr((X,Y)=(0,1))=1pa,Pr((X,Y)=(1,1))=a+p+q1.

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

(1)a=(1p)(1q)+ρpq(1p)(1q).

À 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=q11n(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 0,087n=10p=1/3q=3/4ρ=4/5(1)a=0.003367350.2470.6630.087 ). Voici un tracé de réalisations de la distribution conjointe:1000

Nuage de points

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)
whuber
la source
Cette approche peut-elle être étendue pour générer un nombre quelconque de variables binaires? - pour s'adapter à la matrice de corrélation donnée (ou au plus près pour l'adapter)?
ttnphns
1
2kk2kk1nnn représente.
whuber
C'est un joli résultat. Juste pour revenir un peu sur votre première phrase. Pour obtenir un binôme à partir de variables de Bernoulli indépendantes, n'ont-elles pas besoin d'avoir le même p? Cela n'a aucun effet sur ce que vous avez fait car c'est juste une motivation pour votre approche.
Michael R. Chernick
1
pXqY
@whuber Belle approche! Pouvez-vous me faire savoir s'il y a un document auquel je peux me référer ??
T Nick