2

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

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

4

2 回答 2

1

另一种方法是使用 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())

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

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

于 2020-10-04T20:23:45.027 回答
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)which 计算binom.pmf. 然后减去该数组的第 th 项b和a第 th 项，得到的输出非常接近您的预期。

于 2020-10-03T21:44:46.987 回答