如何从负二项分布生成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)，它从几何分布中生成 10 个随机数，成功概率为 0.5，随机结果是:

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

我测试了它的rNB(3,2,0.1)，它多次从参数为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 和 [=34= 的负二项分布中生成 n 随机数].其中r是成功次数，p是成功概率？

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

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

负二项式distribution:A负二项式实验是一种统计实验，具有以下性质：实验由 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 这个向量。

table 的第一个 r 个数字（通过 table 减去 [1:n] 或等效地 [seq(n)] 得到的是你的负二项式抽奖。我们只是使用 as.numeric 删除 table 的名字。我们还从每个计数中减去成功的次数（即 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)
}

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

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

How to generate n random numbers from negative binomial distribution?

random

r

distribution