python - 从 'partial' 二项分布进行高效采样

Question

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

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

f(x) = bin(x)/sum

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

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

什么是实现这一目标的好/有效的方法？

Answer 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 的累积和。然后将此数组的第 b 项和第 a 项相减，得到的输出非常接近您的预期。