在numpy/scipy中是否有一种函数可以从小对数概率向量中抽取多项式样本，而不会丢失精度？例如：
# sample element randomly from these log probabilities
l = [-900, -1680]
因为下溢，朴素方法会失败：
import scipy
import numpy as np
# this makes a all zeroes
a = np.exp(l) / scipy.misc.logsumexp(l)
r = np.random.multinomial(1, a)
这是一次尝试：
def s(l):
m = np.max(l)
norm = m + np.log(np.sum(np.exp(l - m)))
p = np.exp(l - norm)
return np.where(np.random.multinomial(1, p) == 1)[0][0]
这是最佳/最快的方法吗？最后一步可以避免使用np.exp()吗？
首先，我相信你遇到的问题是因为你错误地对概率进行了归一化处理。以下这行代码是错误的：
a = np.exp(l) / scipy.misc.logsumexp(l)
你正在将一个概率除以一个对数概率，这是没有意义的。相反，你可能想要
a = np.exp(l - scipy.misc.logsumexp(l))
如果您这样做，您会发现a = [1, 0]，并且您的多项式采样器可以按照第二个概率的浮点精度正常工作。
话虽如此，如果您仍需要更高的精度，并且性能不是太大的问题，您可以通过从头开始实现一个多项式采样器，然后修改它以在更高精度上工作。
NumPy的多项式函数是Cython实现的，本质上是对一定数量的二项式样本进行循环，并将它们组合成一个多项式样本。您可以像这样调用它：
np.random.multinomial(10, [0.1, 0.2, 0.7])
# [0, 1, 9]
（注意，这里以及以下的精确输出值都是随机的，并且会在每次调用时发生变化。）
另一种实现多项式采样器的方法是生成N个均匀随机值，然后使用由累积概率定义的直方图来计算它们：
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的增长，它变得非常缓慢。
我思考了一下，意识到有一种更有效的方法，尽管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值上会出现问题——这是因为该算法计算了值(1-p)，所以p的值受到浮点精度的限制。
那么我们该怎么办呢？嗯，事实证明，对于小的p值，Poisson分布是二项式分布的极好近似，并且该实现没有这些问题。因此，我们可以基于一个稳健的二项式采样器构建一个强大的多项式函数，在小的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，多项式分布的近似非常稳健且非常快速。