log(1-exp(x)) 的数值精度

Question

我正在用非常大的数字做一些数学运算（我使用的是 Python，但这个问题不是 Python 特有的）。对于一个值，我有一个公式给我f(t) = Pr(X < t)。我想用这个公式得到Pr(X >= t) = 1 - f(t)。因为f(t) 返回的值非常接近于零，所以我一直在使用对数转换并存储log( f(t) ) 而不是f(t)。我的log( f(t) ) 大约是-1e5 左右。

对于乘法，这非常有效。 log( f(t) * g ) = log( f(t) ) + log(g)。

但是，仅使用 log( f(t) ) 来计算 log( 1 - f(t) ) 非常困难；当然，我可以暂时取我存储和计算 log( 1 - exp( log( f(t) ) ) 的值取幂，但这将返回 log( 1 - 0.0 ) = 0.0，因为 log( f(t) ) 非常接近于零。

你可能会问，“你为什么关心？如果它接近于零，那么 1 减去它就非常接近于 1。”嗯，这是一个很好的观点。你是个聪明的饼干。

问题是我想用它来对值进行排名，所以我真的很关心一个是log(0.999)，另一个是log(0.9999)。您可能还会问：“好吧，为什么不直接对log( f(t) ) 进行排名，然后颠倒顺序以获得log( 1 - f(t) ) 的排名。”再一次，我不得不指出你的问题有多棒。与您交谈真的很愉快。

但问题是：我不只是想按1 - f(t) 排名；我实际上想根据Pr(X >= t) * g(t) = (1 - f(t)) g(t) 进行排名。获取日志后，我得到log( 1 - f(t) ) + log( g(t) )；仅基于f(t) 的排名不会给出正确答案。

过去我写了一个小 Python 函数来计算 log(a + b) 从 log(a) 和 log(b):

def log_add(logA,logB):
    if logA == log(0):
        return logB
    if logA<logB:
        return log_add(logB,logA)
    return log( 1 + math.exp(logB-logA) ) + logA

首先对它们进行归一化以使它们靠近在一起，然后在它们靠近在一起时求幂，这会有所帮助。

不幸的是，我无法使用相同的技巧来进行减法运算，因为没有标准化因子可以将 log(1) 和 log( f(t) ) 靠近在一起，因为它们相距甚远。

有谁知道如何解决这个问题？这似乎是一个经典的问题；我真的希望/希望/祈祷有一个聪明的功能可以在位级别上运行，它可以从log(x) 给我log(1-x)。另外，如果您知道 它是如何工作的，我真的很想知道。

干杯！ 奥利弗

Answer 1

如果log(f(t)) 确实是-1e5（或类似数量级），那么0.0 是log(1-f(t)) 的最佳浮点表示。确实，f(t) = exp(-1e5) 所以，根据 dmuir 提到的泰勒级数，log(1-f(t)) = -exp(-1e5)（这实际上不是一个精确的等式，但它是一个非常好的近似值）。现在，-exp(-1e5) = -3.56e-43430，但是在 0 和 -4e-324 之间没有浮点数，所以最好的浮点表示是 0.0。

因此，使用标准浮点数是不可能做到的。

这有关系吗？你说要根据Pr(X >= t) * g(t) = (1 - f(t)) g(t)排名，相当于log( 1 - f(t) ) + log( g(t) )排名。我们在上面发现了log(1-f(t)) = -3.56e-43430，所以这个术语只有在log(g(t)) 的不同值相差不超过这个小数字并且如果你的计算足够准确以至于它可以通过这些小数字来区分时才会产生影响（如果您使用标准浮点数，那么您的计算将永远不够准确）。换句话说，如果log(f(t)) 确实是-1e5 或类似的值，那么您可以按g(t) 排名。

但是，log(f(t)) 可能是 -1e5 的数量级，但它有时会采用更接近零的值，例如 -10 或 -1。在这种情况下，您不能忽略它，您必须确实按log(1-f(t)) + log(g(t)) 排名。您应该使用math.log1p 函数编写此代码：按log1p(-f(t)) + log(g(t)) 排名。原因是如果 f(t) 接近于零，则 log(1-f(t)) 不准确，但 log1p(-f(t)) 准确。如果 f(t) 非常接近于零，例如当 log(f(t)) = -1e5 时，则为 log1p(-f(t)) = 0.0，因为这是使用标准浮点数所能做到的最好的。

我使用“标准浮点数”是有原因的。可以使用更精确的浮点数，如果您真的想捕获像-3.56e-43430 这样的小数，那您应该这样做。 Python 中的一种可能性是mpmath（不幸的是，它似乎不支持log1p 函数）。请注意，这比标准浮点数要慢得多，正如我所说，我认为您不需要它。但是，如果您想更好地了解这些问题，还是值得一试的。