Python 对数概率求和代码返回天文数字

Question

我有以下代码：

import numpy as np
def sum_log_prob(l):
    # p = l[-1]
    # i = len(l) - 1
    # while i:
    #     i -= 1
    #     curr_p = l[i]
    #     if p > curr_p:
    #         p = p + np.log1p(np.exp(curr_p - p))
    #     else:
    #         p = curr_p + np.log1p(np.exp(p - curr_p))
    p = float('-inf')
    for x in l:
        if p > x:
            p = p + np.log1p(np.exp(x - p))
        else:
            p = x + np.log1p(np.exp(p - x))
        print(p)
    return p

print(sum_log_prob([np.log(0.2), np.log(0.4), np.log(0.2), np.log(0.2)]))

理论上它会取概率的总和，所有这些概率都以对数形式表示。由于对数加法相当于普通乘法，因此我编写了这段代码，以便在概率太小时时最大限度地减少数值下溢。

问题是，它不起作用。这是循环每次迭代的 p：

-1.6094379124341003
-0.5108256237659906
-0.22314355131420965
1.1102230246251565e-16

如您所见，p 非常接近正确答案（0.0，因为 log(1) = log(0.2 + 0.4 + 0.2 + 0.2) = 0），但最后返回一些小数字。

我对这段代码进行了两次实现，都给出了相同的错误答案。其中之一在上面的代码片段中被注释掉了。

请赐教，谢谢！

Answer 1

您面临的问题可能与数值计算中的舍入和浮点精度有关。当您在浮点表示中使用非常小的或非常接近于零的数字时，有可能会丢失精度。

您可以尝试将

np.log1p(np.exp(...))

替换为

np.log(1 + np.exp(...))

以避免出现数值问题。另外，在这种情况下，您可以尝试使用 NumPy 提供的

np.logaddexp

来更安全地计算。

这是更新后的代码：

import numpy as np

def sum_log_prob(l):
    p = float('-inf')
    for x in l:
        p = np.logaddexp(p, x)
    return p

print(sum_log_prob([np.log(0.2), np.log(0.4), np.log(0.2), np.log(0.2)]))

这应该会给出数字上更正确的结果。如果运行此代码，您应该会得到接近于零的结果，正如您所期望的那样。