'partial' 二项分布的高效抽样

Question

我想从二项分布 B(n,p) 中采样，但有一个附加约束，即采样值属于范围 [a,b]（而不是正常的 0 到 n 范围）。换句话说，我必须从二项分布中抽取一个值，因为它位于 [a,b] 范围内。在数学上，我可以根据二项分布 bin(x) = [(nCx)*(p)^x*(1-p)^(n-x)] 的 pmf 将此分布 (f(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 然后使用统一随机变量生成样本的方法似乎效率极低.

实现此目标的 nice/efficient 方法是什么？

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 项给出的输出非常接近您的预期。

Answer 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())

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

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

'partial' 二项分布的高效抽样

Efficient sampling from a 'partial' binomial distribution

python

random

numpy

binomial-cdf

错误修复