python - 从小对数概率向量中以numpy/scipy采样多项式

Question

numpy / scipy中是否有一个函数可以让您从对数概率小的 vector 中采样多项式而又不损失精度？例:

# 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()吗？

Answer 1

首先，我相信您遇到的问题是因为您不正确地归一化了概率。这行是不正确的:

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!请注意，我们不能使用大于 N的 1E18，因为它将溢出长整数。
但是我们可以使用 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的多项式分布的非常鲁棒且非常快速的近似值。