Menghasilkan variabel acak binomial berkorelasi

Saya bertanya-tanya apakah mungkin untuk menghasilkan variabel binomial acak berkorelasi mengikuti pendekatan transformasi linier?

Di bawah ini, saya mencoba sesuatu yang sederhana dalam R dan menghasilkan beberapa korelasi. Tapi saya bertanya-tanya apakah ada cara berprinsip untuk melakukan ini?

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)

Variabel binomial biasanya dibuat dengan menjumlahkan variabel Bernoulli independen. Mari kita lihat apakah kita bisa mulai dengan sepasang variabel Bernoulli yang berkorelasi dan melakukan hal yang sama.$(X,Y)$

Misalkan adalah variabel Bernoulli ( p ) (yaitu, Pr ( X = 1 ) = p dan$X$$(p)$$\Pr(X=1)=p$ ) dan Y adalahvariabelBernoulli ( q ) . Untuk menjabarkan distribusi bersama mereka, kita perlu menentukan keempat kombinasi hasil. Penulisan Pr ( ( X , Y ) = ( 0 , 0 ) ) $\Pr(X=0)=1-p$$Y$$(q)$ kita dapat dengan mudah mengetahui sisanya dari aksioma probabilitas: Pr ( ( X , Y ) = ( 1 , 0 ) ) = 1 - q - a

Memasukkan ini ke dalam rumus untuk koefisien korelasi dan penyelesaiannya memberikan a = ( 1 - p ) ( 1 - q ) + ρ $\rho$

Asalkan keempat probabilitas non-negatif, ini akan memberikan distribusi gabungan yang valid - dan solusi ini membuat parameter semua distribusi Bernoulli bivariat. (Ketika , ada solusi untuk semua korelasi yang bermakna secara matematis antara - 1 dan 1. ) Ketika kita menjumlahkan n dari variabel-variabel ini, korelasinya tetap sama - tetapi sekarang distribusi marjinal adalah Binomial ( n , p ) dan Binomial ( n , q ) , seperti yang diinginkan.$p=q$$-1$$1$$n$$(n,p)$$(n,q)$

Contoh

Biarkan , p = 1 / 3 , q = 3 / 4 , dan kami ingin korelasi menjadi ρ = - $n=10$$p=1/3$$q=3/4$$\rho=-4/5$$(1)$$a=0.00336735$$0.247$$0.663$$0.087$$1000$

Garis merah menunjukkan rata-rata sampel dan garis putus-putus adalah garis regresi. Mereka semua dekat dengan nilai yang dimaksudkan. Poin-poin telah dikelompokkan secara acak dalam gambar ini untuk menyelesaikan tumpang tindih: setelah semua, distribusi Binomial hanya menghasilkan nilai-nilai integral, sehingga akan ada banyak overplotting.

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

Kode

Berikut ini adalah Rimplementasinya.

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