从numpy/scipy中的小对数概率向量中抽取多项式

Question

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()在最后一步避免吗？

Answer 1

首先,我相信您遇到的问题是因为您错误地将您的概率标准化.这行不正确:

a = np.exp(l) / scipy.misc.logsumexp(l)

您将概率除以对数概率,这没有任何意义.相反,你可能想要

a = np.exp(l - scipy.misc.logsumexp(l))

如果这样做,您会发现a = [1, 0]并且您的多项式采样器按预期工作,在第二个概率中达到浮点精度.

小N:直方图的解决方案

也就是说,如果你仍然需要更高的精度和性能并不是一个问题,你可以取得进步的一种方法是从头开始实现多项式采样器,然后修改它以更高的精度工作.

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:泊松近似的一种解法

上述解决方案的问题在于,随着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值,泊松分布是二项分布的非常好的近似,并且实现没有这些问题.因此,我们可以构建一个强大的多项式函数,它基于一个强大的二项式采样器,可以在小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.