如何从负二项分布中生成 n 个随机数？

Question

我正在尝试创建一个函数，以便从负二项分布中生成 n 随机数。 为了生成它，我首先制作了一个函数来从几何分布中生成n 随机变量。我从几何分布生成n随机数的函数如下：

rGE<-function(n,p){
  I<-rep(NA,n)
  for (j in 1:n){
  x<-rBer(1,p)
  i<-1 # number of trials
  while(x==0){
    x<-rBer(1,p)
    i<-i+1
  }
  I[j]<- i
  }
  return(I)
}

我测试了这个函数（rGE），例如rGE(10,0.5)，它从成功概率为0.5的几何分布中生成10随机数，随机结果是：

[1] 2 4 2 1 1 3 4 2 3 3

在rGE 函数中，我使用了一个名为rBer 的函数，即：

rBer<-function(n,p){
  sample(0:1,n,replace = TRUE,prob=c(1-p,p))
}

现在，我想改进上面的函数 (rGE)，以便创建一个函数来从负二项式函数生成 n 随机数。我做了以下功能：

rNB<-function(n,r,p){
  I<-seq(n)
  for (j in 1:n){
    x<-0
    x<-rBer(1,p)
    i<-1 # number of trials
    while(x==0 & I[j]!=r){
      x<-rBer(1,p)
      i<-i+1
    }
    I[j]<- i
  }
  return(I)
}

我对@9​​87654335@ 进行了多次测试，它从带有参数r=2 和p=0.1 的负二项分布中生成了3 个随机数：

> rNB(3,2,0.1)
[1] 2 1 7
> rNB(3,2,0.1)
[1] 3 1 4
> rNB(3,2,0.1)
[1] 3 1 2
> rNB(3,2,0.1)
[1] 3 1 3
> rNB(3,2,0.1)
[1] 46  1 13 

如您所见，我认为我的函数 (rNB) 无法正常工作，因为结果总是为第二个随机数生成 1。 谁能帮我纠正我的函数（rNB），以便从参数为n、r 和p 的负二项分布中生成n 随机数。其中r是成功次数，p是成功概率？

[[提示：关于几何分布和负二项分布的解释： 几何分布：在概率论和统计学中，几何分布是两种离散概率分布中的一种：

1. 获得一次成功所需的伯努利试验次数 X 的概率分布，支持集合 {1,2,3,...}。
2. 第一次成功前失败次数 Y = X − 1 的概率分布，支持集合 { 0, 1, 2, 3, ... }

负二项分布：负二项实验是一种统计实验，具有以下属性： 该实验由 x 次重复试验组成。 每个试验只能产生两种可能的结果。我们称其中一个结果为成功，另一个为失败。 以 P 表示的成功概率在每次试验中都是相同的。 试验是独立的；也就是说，一项试验的结果不会影响其他试验的结果。 实验继续进行，直到观察到 r 次成功，其中 r 是预先指定的。 ]]

Answer 1

如果你使用 R 的原生向量化，你的函数会更快。这样做的方法是一次生成所有伯努利试验。

请注意，对于负二项分布，预期值（即获得r 成功所需的平均伯努利试验次数）为r * p / (1 - p) (Reference)

如果我们要抽取n 负二项式样本，那么伯努利试验的预期总数将因此为n * r * p / (1 - p)。所以我们想至少抽取那么多伯努利样本。为简单起见，我们可以从绘制两倍的数字开始：2 * n * r * p / (1 - p)。万一这还不够，我们可以再次重复绘制两倍，直到我们有足够的；一旦伯努利试验的合成向量之和大于r * n，我们就知道我们有足够的伯努利试验来模拟我们的n 负二项试验。

我们现在可以在伯努利试验的向量上运行cumsum 来跟踪阳性试验的数量。如果您随后通过%/% r 对该向量执行整数除法，您将根据它们所属的负二项式试验标记所有伯努利试验。然后你table这个向量。

表格的第一个r 数字（通过将表格子集[1:n] 或等效地通过[seq(n)] 获得是您的负二项式抽签。我们只是使用as.numeric 删除表格的名称。我们还减去成功的数量（即r），来自我们的每个计数，因为我们只计算失败，而不是成功。

rNB <- function(n, r, p) {
  mult <- 2
  all_samples <- 0
  while(sum(all_samples) < n * r)
  {
    all_samples <- rBer(mult * n * r * p / (1 - p), p)
    mult <- mult * 2
  }
  as.numeric(table(cumsum(all_samples) %/% r))[seq(n)] - r
}

所以我们可以这样做：

rNB(3, 2, 0.1)
#> [1] 14 19 41

rNB(3, 2, 0.1)
#> [1] 23  6 56

rNB(3, 2, 0.1)
#> [1] 11 31 59

rNB(3, 2, 0.1)
#> [1]  7 21 14

mean(rNB(10000, 2, 0.1))
#> [1] 18.0002

我们可以针对 R 自己的 rnbinom 进行测试：

mean(rnbinom(10000, 2, 0.1))
#> [1] 18.0919

hist(rnbinom(10000, 2, 0.5), breaks = 0:20)

hist(rNB(10000, 2, 0.5), breaks = 0:20)

请注意，您自己版本的逻辑并不完全正确。特别是，while(x == 0 & I[j] != r) 行没有任何意义。 I 是 1:n 的向量，因此在您的示例中，每当 j 为 2 时，I[j] 等于 r 并且循环停止。这就是为什么你的第二个数字总是 1。我不知道你在这里想做什么。

如果您想一次进行一次伯努利试验，就像您在自己的版本中所做的那样，请尝试这个修改后的功能。变量名应该可以让逻辑更容易理解：

rNB <- function(n, r, p) {
  # Create an empty vector of length n for our results
  draws <- numeric(n)
  
  # Now for each of the n trials we will get a negative binomial sample:
  for (i in 1:n) {
    # Create success and failure counters for this draw
    failures  <- successes <- 0
    
    # Now run Bernoulli trials, counting successes and failures as we go
    # until we hit r successes
    while(successes < r)
    {
      if(rBer(1, p) == 1) 
        successes <- successes + 1
      else
        failures  <- failures + 1
    }

    # Once we have reached r successes, the current number of failures is our
    # negative binomial draw
    draws[i] <- failures
  }
  
  return(draws)
}

这与更快但更不透明的矢量化版本给出了相同的结果。