如何在给定一系列 bin 概率的情况下生成 bin 计数的随机样本？答案

【问题标题】：How can I generate a random sample of bin counts given a sequence of bin probabilities?如何在给定一系列 bin 概率的情况下生成 bin 计数的随机样本？
【发布时间】：2015-10-22 04:19:11
【问题描述】：

我有一个整数，需要根据概率分布分成多个 bin。例如，如果我有N=100 对象进入[0.02, 0.08, 0.16, 0.29, 0.45]，那么你可能会得到[1, 10, 20, 25, 44]。

import numpy as np
# sample distribution
d = np.array([x ** 2 for x in range(1,6)], dtype=float)
d = d / d.sum()
dcs = d.cumsum()
bins = np.zeros(d.shape)
N = 100
for roll in np.random.rand(N):
    # grab the first index that the roll satisfies
    i = np.where(roll < dcs)[0][0]  
    bins[i] += 1

实际上，N 和我的 bin 数量非常大，因此循环并不是一个真正可行的选择。有什么方法可以矢量化此操作以加快速度？

【问题讨论】：

标签： python numpy random vectorization probability-density

【解决方案1】：

您可以通过获取 cumsum 将 PDF 转换为 CDF，使用它来定义一组介于 0 和 1 之间的 bin，然后使用这些 bin 计算 N 长随机数的直方图统一向量：

cdf = np.cumsum([0, 0.02, 0.08, 0.16, 0.29, 0.45])     # leftmost bin edge = 0
counts, edges = np.histogram(np.random.rand(100), bins=cdf)

print(counts)
# [ 4,  8, 16, 30, 42]

【讨论】：

当然，我使用术语 binning 并忽略了查找直方图。非常感谢！
不用担心。我冒昧地编辑了您问题的标题，以便其他人更容易找到。

【解决方案2】：

您可以使用np.bincount 进行分箱操作，同时使用np.searchsorted 来执行等效于roll < dcs 的操作。这是实现这些承诺的实现 -

bins = np.bincount(np.searchsorted(dcs,np.random.rand(N),'right'))

使用给定参数的运行时测试 -

In [72]: %%timeit
    ...: for roll in np.random.rand(N):
    ...:     # grab the first index that the roll satisfies
    ...:     i = np.where(roll < dcs)[0][0]  
    ...:     bins[i] += 1
    ...: 
1000 loops, best of 3: 721 µs per loop

In [73]: %%timeit
    ...: np.bincount(np.searchsorted(dcs,np.random.rand(N),'right'))
    ...: 
100000 loops, best of 3: 13.5 µs per loop

【讨论】：

对于较小的系统大小，您的解决方案比其他解决方案快 10 倍，但对于较大的系统大小，速度会慢 10 倍。很好，感谢您提供替代解决方案！

【解决方案3】：

另一种方法：

import numpy as np

p = [0.02, 0.08, 0.16, 0.29, 0.45]
np.bincount(np.random.choice(range(len(p)), size=100, p=p), minlength=len(p))
# array([ 1,  6, 16, 25, 52])

似乎没有必要分配一个长度为 100 的数组，但我还没有在 numpy 中看到避免它的方法。

【讨论】：