首先，我相信您遇到的问题是因为您错误地将概率标准化。此行不正确：
a = np.exp(l) / scipy.misc.logsumexp(l)
您将概率除以对数概率，这是没有意义的。相反，您可能想要
a = np.exp(l - scipy.misc.logsumexp(l))
如果您这样做，您会找到a = [1, 0]，并且您的多项式采样器可以按预期工作，直到第二个概率的浮点精度。
小 N 的解决方案：直方图
也就是说，如果您仍然需要更高的精度并且性能不是那么重要，那么您可以取得进展的一种方法是从头开始实施多项式采样器，然后对其进行修改以提高精度。
NumPy 的多项式函数是implemented in Cython，本质上是对多个二项式样本执行循环，并将它们组合成一个多项式样本。
你可以这样称呼它：
np.random.multinomial(10, [0.1, 0.2, 0.7])
# [0, 1, 9]
（请注意，此处和以下的精确输出值是随机的，并且会因调用而异）。
另一种实现多项式采样器的方法是生成 N 个均匀随机值，然后使用由累积概率定义的 bin 计算直方图：
def multinomial(N, p):
rand = np.random.uniform(size=N)
p_cuml = np.cumsum(np.hstack([[0], p]))
p_cuml /= p_cuml[-1]
return np.histogram(rand, bins=p_cuml)[0]
multinomial(10, [0.1, 0.2, 0.7])
# [1, 1, 8]
考虑到这种方法，我们可以考虑通过将所有内容保存在日志空间中来实现更高的精度。主要技巧是要意识到均匀随机偏差的对数等效于指数随机偏差的负数，因此您可以在不留对数空间的情况下完成上述所有操作：
def multinomial_log(N, logp):
log_rand = -np.random.exponential(size=N)
logp_cuml = np.logaddexp.accumulate(np.hstack([[-np.inf], logp]))
logp_cuml -= logp_cuml[-1]
return np.histogram(log_rand, bins=logp_cuml)[0]
multinomial_log(10, np.log([0.1, 0.2, 0.7]))
# [1, 2, 7]
即使对于 p 数组中非常小的值，生成的多项式绘制也将保持精度。
不幸的是，这些基于直方图的解决方案将比原生 numpy.multinomial 函数慢很多，因此如果性能是一个问题，您可能需要另一种方法。一种选择是将上面链接的 Cython 代码调整为在对数空间中工作，使用与我在这里使用的类似的数学技巧。
大 N 的解决方案：泊松近似
上述解决方案的问题是随着N变大，它变得非常慢。
我正在考虑这一点，并意识到有一种更有效的前进方式，尽管np.random.multinomial 失败的概率小于1E-16 左右。
这是一个失败的例子：在 64 位机器上，由于代码的实现方式，第一个条目总是为零，而实际上它应该给出接近 10 的值：
np.random.multinomial(1E18, [1E-17, 1])
# array([ 0, 1000000000000000000])
如果您深入研究源代码，您可以将此问题追溯到构建多项式函数的二项式函数。 cython 代码在内部执行如下操作：
def multinomial_basic(N, p, size=None):
results = np.array([np.random.binomial(N, pi, size) for pi in p])
results[-1] = int(N) - results[:-1].sum(0)
return np.rollaxis(results, 0, results.ndim)
multinomial_basic(1E18, [1E-17, 1])
# array([ 0, 1000000000000000000])
问题在于binomial 函数在p 的非常小的值上阻塞——这是因为computes the value (1 - p) 算法，所以p 的值受到浮点精度的限制。
那么我们能做些什么呢？好吧，事实证明，对于较小的 p 值，泊松分布是二项分布的极好的近似值，并且实现不存在这些问题。因此，我们可以基于一个稳健的二项式采样器构建一个稳健的多项式函数，该采样器在小 p 处切换到泊松采样器：
def binomial_robust(N, p, size=None):
if p < 1E-7:
return np.random.poisson(N * p, size)
else:
return np.random.binomial(N, p, size)
def multinomial_robust(N, p, size=None):
results = np.array([binomial_robust(N, pi, size) for pi in p])
results[-1] = int(N) - results[:-1].sum(0)
return np.rollaxis(results, 0, results.ndim)
multinomial_robust(1E18, [1E-17, 1])
array([ 12, 999999999999999988])
第一个条目是非零且接近 10，正如预期的那样！请注意，我们不能使用大于1E18 的N，因为它会溢出长整数。
但是我们可以使用size 参数确认我们的方法适用于较小的概率，并对结果进行平均：
p = [1E-23, 1E-22, 1E-21, 1E-20, 1]
size = int(1E6)
multinomial_robust(1E18, p, size).mean(0)
# array([ 1.70000000e-05, 9.00000000e-05, 9.76000000e-04,
# 1.00620000e-02, 1.00000000e+18])
我们看到，即使对于这些非常小的概率，多项式值也以正确的比例出现。结果是对小p 的多项分布的非常稳健且非常快速的近似。