从R中的二项分布生成相关随机数

Question

我试图找到一种从几个二项分布生成相关随机数的方法.

我知道如何使用正态分布(使用mvrnorm),但我找不到适用于二项式的函数.

Answer 1

您可以使用copula包生成相关的制服,然后使用该qbinom函数将它们转换为二项式变量.这是一个简单的例子:

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

现在x2是一个矩阵,其中2列代表2个相关的二项式变量.

Answer 2

在每次试验中具有n次试验和成功概率p的二项式变量可以被视为n个伯努利试验的总和,每个试验也具有成功概率p.

类似地,您可以通过对具有所需相关性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)

值得注意的是,有许多不同的关节分布共享所需的相关系数.模拟方法rmvBinomial()生成其中一个,但它是否合适将取决于生成数据的过程.

正如这个R-help对类似问题的回答所说的那样(然后继续更详细地解释这个想法):

虽然二元正态(给定平均值和方差)由相关系数唯一定义,但二元二项式不是这种情况.