从“部分”二项分布中进行有效抽样答案

【问题标题】：Efficient sampling from a 'partial' binomial distribution从“部分”二项分布中进行有效抽样
【发布时间】：2020-10-03 19:23:08
【问题描述】：

我想从二项分布 B(n,p) 中采样，但有一个额外的约束，即采样值属于 [a,b] 范围（而不是正常的 0 到 n 范围）。换句话说，我必须从二项分布中采样一个值，因为它位于 [a,b] 范围内。在数学上，我可以将这个分布的 pmf (f(x)) 写成二项分布 bin(x) = [(nCx)*(p)^x*(1-p)^(n-x)] 的 pmf 为

sum = 0
for i in range(a,b+1):
    sum += bin(i)

f(x) = bin(x)/sum

从该分布中采样的一种方法是采样一个均匀分布的数字并应用 CDF 的倒数（使用 pmf 获得）。但是，我认为这不是一个好主意，因为 pmf 计算很容易变得非常耗时。

在我的情况下n,x,a,b 的值非常大，由于nCx 中的阶乘项，这种计算 pmf 然后使用统一随机变量生成样本的方式似乎效率极低。

有什么好的/有效的方法来实现这一点？

【问题讨论】：

我认为二项分布的范围必须是从 0 到 n。你能为我提供这个问题背后的进一步数学的链接吗？我想了解更多。
另外，如果您在优化代码时遇到问题，我建议您使用 NumPy 来向量化您的计算。它比 for 循环快 100 倍
我找到了一个无错误且有效的解决方案，它允许为数百万个项目计算 pmf 和 cdf。它需要使用scipy.stats.binom。这对我来说是新事物，因为到目前为止我只知道 scipy.special.comb 具有误导性。
@KevinChoi 我在我的项目中模拟二项式 RV x 以获得另一个 RV y，它是 x 的函数。但是，y 也有范围 [0,n]。为了确保y = f(x)在可接受的范围内，需要对x施加一个条件，归结为在[a,b]之间二项式选择x。

标签： python numpy random binomial-cdf

【解决方案1】：

这是一种在很短的时间内收集bin 的所有值的方法：

from scipy.special import comb
import numpy as np
def distribution(n, p=0.5):
    x = np.arange(n+1)
    return comb(n, x, exact=False) * p ** x * (1 - p) ** (n - x)

n=1000 可以在四分之一微秒内完成。

示例运行：

>>> distribution(4):
array([0.0625, 0.25  , 0.375 , 0.25  , 0.0625])

您可以像这样对该数组的特定部分求和：

>>> np.sum(distribution(4)[2:4])
0.625

备注：对于n>1000，这个分布的中间值需要在乘法中使用非常大的数字，因此RuntimeWarning被提高了。

错误修正

您可以等效地使用scipy.stats.binom：

from scipy.stats import binom
def distribution(n, p):
    return binom.pmf(np.arange(n+1), n, p)

这与上述方法非常有效（n=1000000 在三分之一秒内）。或者，您可以使用binom.cdf(np.arange(n+1), n, p) 计算binom.pmf 的累积总和。然后将该数组的 bth 和 ath 项相减，得到的输出非常接近您的预期。

【讨论】：

哇，这似乎真的很快。我想知道他们如何在binom 类中如此高效地计算nCr。

【解决方案2】：

另一种方法是使用 CDF，它是逆向的，例如：

from scipy import stats

dist = stats.binom(100, 0.5)

# limit ourselves to [60, 100]
lo, hi = dist.cdf([60, 100])

# draw a sample
x = dist.ppf(stats.uniform(lo, hi-lo).rvs())

应该给我们范围内的值。请注意，由于浮点精度，这可能会给您提供超出您想要的值的值。它在分布的平均值之上变得更糟

请注意，对于较大的值，您不妨使用正态近似值

【讨论】：