来自 R 中的自定义分布的样本答案

【问题标题】：Sample from custom distribution in R来自 R 中的自定义分布的样本
【发布时间】：2017-08-13 22:39:32
【问题描述】：

我已经在 R 中实现了负二项分布的替代参数化，就像这样（另见 here）：

nb = function(n, l, a){
  first = choose((n + a - 1), a-1)
  second = (l/(l+a))^n
  third = (a/(l+a))^a
  return(first*second*third)
}

其中 n 是计数，lambda 是平均值，a 是过度离散项。

我想从此分布中抽取随机样本，以验证我对负二项式混合模型的实现，但不知道如何去做。这个函数的 CDF 不容易定义，所以我考虑尝试像 here 讨论的那样尝试拒绝抽样，但这也不起作用（我不确定为什么 - 文章说首先从两者之间的均匀分布中提取0 和 1，但我希望我的 NB 分布对整数计数进行建模...我不确定我是否完全理解这种方法。）

感谢您的帮助。

【问题讨论】：

标签： r statistics distribution sampling

【解决方案1】：

看来你可以：

1) 画一个介于 0 和 1 之间的均匀随机数。

2) 对概率密度函数进行数值积分（这实际上只是一个和，因为分布是离散的且下界为零）。

3) 无论您的集成中的哪个值使 cdf 超过您的随机数，这就是您的随机抽奖。

所以大家一起做以下事情：

r <- runif(1,0,1)
cdf <- 0
i <- -1
while(cdf < r){
  i <- i+1
  p <- PMF(i)
  cdf <- cdf + p
}

其中 PMF(i) 是 i 计数上的概率质量，由分布参数指定。此 while 循环结束时 i 的值就是您的示例。

【讨论】：

这似乎类似于我在帖子中链接到的拒绝抽样方法（也许我弄错了）？但我会试一试。
这与拒绝抽样有很大不同。拒绝抽样的工作原理是提出一个样本，计算似然性，并以与似然成比例的概率接受样本。（从 0 到 1 的均匀分布出现在链接拒绝抽样示例中，因为 beta 分布是在该区间上定义的，因此这就是建议样本来自的区间）。我的答案是找到与 cdf 的给定值对应的样本，其中给定值均匀分布在 0 和 1（cdf 的端点）之间。
这里的重点是，使用这个停止标准进行积分（求和）会计算逆 cdf。而且因为分布是离散的，所以计算是精确的。这不是一个封闭的形式，所以如果你需要大量的抽奖，你会付出效率的代价。但它仍然应该很快，尤其是当您从中采样的分布的预期值相对较小时。
Jacob，感谢您的详细解答！我使用了你的方法，得到了一个看起来相当合理的采样密度。

【解决方案2】：

我建议您查看 Uniform 分布以及 Uniform 的普遍性。您可以通过将一个均匀分布的变量传递给 NB 二项式的逆 CDF 来做您想做的事情，您将得到的是从您的 NB 二项式分布中采样的一组点。

编辑：我看到负二项式的 CDF 没有封闭形式的逆。我的第二个建议是废弃您的函数并使用内置函数：

library(MASS)
rnegbin(n, mu = n, theta = stop("'theta' must be specified"))

【讨论】：

感谢您的评论，但我想保留我的功能，因为 NB 的这种公式通常用于基因表达/读取计数，这是我正在尝试建模的系统。
好吧，我建议您开始为您的表示搜索一个封闭形式的表达式，然后使用统一的普遍性。
在这种情况下，您实际上并不需要一个封闭的表格，除非效率是一个主要问题。用数值计算逆 CDF 很容易。而且由于分布是离散的且有界为零，因此计算是精确的。

【解决方案3】：

如果您真的只是想测试并且速度不是问题，那么正如其他人提到的，反转方法可能是要走的路。

对于离散随机变量，它需要一个简单的 while 循环。看 Non-Uniform Random Variate Generation L. Devroye，第 3 章，p. 85.

【讨论】：