Solución:
Puede generar uniformes correlacionados utilizando el copula paquete, luego use el qbinom función para convertirlos en variables binomiales. Aquí hay un ejemplo rápido:
library(copula)
tmp <- normalCopula( 0.75, dim=2 )
x <- rcopula(tmp, 1000)
x2 <- cbind( qbinom(x[,1], 10, 0.5), qbinom(x[,2], 15, 0.7) )
Ahora x2 es una matriz con 2 columnas que representan 2 variables binomiales que están correlacionadas.
Una variable binomial con n ensayos y probabilidad p de éxito en cada ensayo puede verse como la suma de n ensayos de Bernoulli, cada uno de los cuales también tiene probabilidad p de éxito.
De manera similar, puede construir pares de variables binomiales correlacionadas sumando pares de variables de Bernoulli que tengan la correlación deseada r.
require(bindata)
# Parameters of joint distribution
size <- 20
p1 <- 0.5
p2 <- 0.3
rho<- 0.2
# Create one pair of correlated binomial values
trials <- rmvbin(size, c(p1,p2), bincorr=(1-rho)*diag(2)+rho)
colSums(trials)
# A function to create n correlated pairs
rmvBinomial <- function(n, size, p1, p2, rho) {
X <- replicate(n, {
colSums(rmvbin(size, c(p1,p2), bincorr=(1-rho)*diag(2)+rho))
})
t(X)
}
# Try it out, creating 1000 pairs
X <- rmvBinomial(1000, size=size, p1=p1, p2=p2, rho=rho)
# cor(X[,1], X[,2])
# [1] 0.1935928 # (In ~8 trials, sample correlations ranged between 0.15 & 0.25)
Es importante tener en cuenta que hay muchas distribuciones conjuntas diferentes que comparten el coeficiente de correlación deseado. El método de simulación en rmvBinomial() produce uno de ellos, pero si es el apropiado o no dependerá del proceso que esté generando sus datos.
Como se indica en esta respuesta de R-help a una pregunta similar (que luego explica la idea con más detalle):
mientras que una normal bivariada (medias y varianzas dadas) se define de forma única por el coeficiente de correlación, este no es el caso de un binomio bivariado