Gerando Variáveis ​​Aleatórias Binomiais Correlacionadas

21

Eu queria saber se seria possível gerar variáveis ​​binomiais aleatórias correlacionadas seguindo uma abordagem de transformação linear?

Abaixo, tentei algo simples em R e isso produz alguma correlação. Mas eu queria saber se existe uma maneira de fazer isso com princípios?

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
fonte
2
Y 2 X 1 = X 2 = 1 Y 1 = 1,5 Y iY1 e podem estar correlacionados, mas não serão mais binomiais. Exemplo, e portanto, o não pode ser variável aleatória binomial. Eu sugiro que você analise a distribuição Multinomial. Y2X1=X2=1Y1=1.5Yi
knrumsey - Restabelece Monica
1
A resposta curta para a pergunta é procurar a palavra-chave copula, que ajuda na geração de variáveis ​​dependentes com margens fixas.
Xian

Respostas:

32

As variáveis ​​binomiais são geralmente criadas pela soma de variáveis ​​independentes de Bernoulli. Vamos ver se podemos começar com um par de variáveis ​​Bernoulli correlacionadas e fazer a mesma coisa.(X,Y)

Suponhamos que é um Bernoulli variável (isto é, e ) e é um Bernoulli variável. Para determinar sua distribuição conjunta, precisamos especificar todas as quatro combinações de resultados. Escrevendo podemos facilmente descobrir o resto a partir dos axiomas da probabilidade:( p ) Pr ( X = 1 ) = pX(p)Pr(X=1)=pY ( q ) Pr ( ( X , Y ) = ( 0 , 0 ) ) = a , Pr ( ( X , Y ) = ( 1 , 0 ) ) = 1 - q - a ,Pr(X=0)=1pY(q)

Pr((X,Y)=(0,0))=a,
Pr((X,Y)=(1,0))=1qa,Pr((X,Y)=(0,1))=1pa,Pr((X,Y)=(1,1))=a+p+q1.

Ao inserir isso na fórmula do coeficiente de correlação e resolver, obtemosρ

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

Desde que todas as quatro probabilidades não sejam negativas, isso fornecerá uma distribuição conjunta válida - e essa solução parametriza todas as distribuições bivariadas de Bernoulli. (Quando , existe uma solução para todas as correlações matematicamente significativas entre e ) Quando somamos dessas variáveis, a correlação permanece a mesma - mas agora as distribuições marginais são Binomial e Binomial , conforme desejado.- 1 1 np=q11n( n , q )(n,p)(n,q)

Exemplo

Seja , , , e gostaríamos que a correlação fosse . A solução para é (e as outras probabilidades estão em torno de , e ). Aqui está um gráfico de realizações da distribuição conjunta:p = 1 / 3 q = 3 / 4 ρ = - 4 / 5 ( 1 ) um = 0,00336735n=10p=1/3q=3/4ρ=4/5(1)a=0.003367350.2470.6630.0871000

Gráfico de dispersão

As linhas vermelhas indicam as médias da amostra e a linha pontilhada é a linha de regressão. Todos estão próximos dos valores pretendidos. Os pontos foram aleatoriamente alterados nesta imagem para resolver as sobreposições: afinal, as distribuições binomiais produzem apenas valores integrais, portanto haverá uma grande quantidade de plotagem.

Uma maneira de gerar essas variáveis ​​é amostrar vezes de com as probabilidades escolhidas e depois converter cada em , cada em , cada em e cada em . Soma os resultados (como vetores) para obter uma realização 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

Aqui está uma Rimplementação.

#
# 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
fonte
Essa abordagem pode ser estendida para gerar qualquer número de variáveis ​​binárias? - ajustar-se à matriz de correlação fornecida (ou no máximo próximo a ajustá-la)?
ttnphns
1
@ttnphns Sim, mas as opções explodem: a tabela de probabilidades deve ser determinada por parâmetros marginais, a restrição de soma para unidade e (portanto) parâmetros adicionais. Evidentemente, você teria muita liberdade para selecionar (ou restringir) esses parâmetros, de acordo com as propriedades multivariadas que deseja criar. Além disso, uma abordagem semelhante pode ser usada para gerar variáveis ​​binomiais correlacionadas com diferentes valores de seus parâmetros " ". Parvin: Acredito que "quando somamos dessas variáveis" explica sem ambiguidade o que representa. k 2 k - k - 1 n n n2kk2kk1nnn
whuber
Este é um bom resultado. Só para escolher um pouco da sua primeira frase. Para obter um binômio a partir de variáveis ​​independentes de Bernoulli, eles não precisam ter o mesmo p? Isso não afeta o que você fez, pois é apenas uma motivação para você se aproximar.
Michael R. Chernick
1
@ Michael Obrigado - você está certo. Outro motivo pelo qual não se relaciona com o método descrito aqui é que esse método ainda envolve somar variáveis ​​de Bernoulli com um parâmetro comum: o parâmetro é para todas as variáveis e para todas as variáveisPara manter o post razoavelmente curto, não apresentei histogramas das distribuições marginais para demonstrar que eles são realmente binomiais (mas na verdade fiz isso na minha análise original apenas para garantir que eles estavam funcionando!). X q YpXqY
whuber
@whuber Boa abordagem! Você pode me informar se há algum artigo ao qual possa me referir?
T Nick