Generando variables aleatorias binomiales correlacionadas

Me preguntaba si sería posible generar variables binomiales aleatorias correlacionadas siguiendo un enfoque de transformación lineal.

A continuación, probé algo simple en R y produce cierta correlación. Pero me preguntaba si hay una forma de principios para hacer esto.

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)

Las variables binomiales generalmente se crean sumando variables independientes de Bernoulli. Veamos si podemos comenzar con un par de variables de Bernoulli correlacionadas y hacer lo mismo.$(X,Y)$

Supongamos que es una variable de Bernoulli (es decir, y ) e es una variable de Bernoulli . Para precisar su distribución conjunta, necesitamos especificar las cuatro combinaciones de resultados. Escribiendo podemos calcular fácilmente el resto a partir de los axiomas de probabilidad:$X$$(p)$Pr ( X = 0 ) = 1 - p Y ( q ) Pr ( ( X , Y ) = ( 0 , 0 ) ) = a , Pr ( ( X , Y ) = ( 1 , 0 ) ) = 1 - $\Pr(X=1)=p$$\Pr(X=0)=1-p$$Y$$(q)$

Al conectar esto a la fórmula para el coeficiente de correlación y resolver, se obtienea = ( 1 - p ) ( 1 - q ) + ρ $\rho$

Siempre que las cuatro probabilidades no sean negativas, esto proporcionará una distribución conjunta válida, y esta solución parametriza todas las distribuciones bivariadas de Bernoulli. (Cuando , hay una solución para todas las correlaciones matemáticamente significativas entre y ) Cuando sumamos de estas variables, la correlación sigue siendo la misma, pero ahora las distribuciones marginales son Binomiales y Binomial , según se desee.- 1 1 n ( n , p ) ( n , q $p=q$$-1$$1$$n$$(n,p)$$(n,q)$

Ejemplo

Sea , , , y nos gustaría que la correlación sea . La solución a es (y las otras probabilidades son alrededor de , y ). Aquí hay una gráfica de realizaciones de la distribución conjunta:p = 1 / 3 q = 3 / 4 ρ = - 4 / 5 ( 1 ) un = 0,00336735 0,247 0,663 0,087 $n=10$$p=1/3$$q=3/4$$\rho=-4/5$$(1)$$a=0.00336735$$0.247$$0.663$$0.087$$1000$

Las líneas rojas indican las medias de la muestra y la línea de puntos es la línea de regresión. Todos están cerca de sus valores previstos. Los puntos se han alterado aleatoriamente en esta imagen para resolver las superposiciones: después de todo, las distribuciones binomiales solo producen valores integrales, por lo que habrá una gran cantidad de superposición.

Una forma de generar estas variables es muestrear veces desde con las probabilidades elegidas y luego convertir cada en , cada en , cada en , y cada en . Suma los resultados (como vectores) para obtener una realización de .{ 1 , 2 , 3 , 4 } 1 ( 0 , 0 ) 2 ( 1 , 0 ) 3 ( 0 , 1 ) 4 ( 1 , 1 ) ( X , Y $n$$\{1,2,3,4\}$$1$$(0,0)$$2$$(1,0)$$3$$(0,1)$$4$$(1,1)$$(X,Y)$

Código

Aquí hay una Rimplementación.

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