在python中，如何计算1减去给定小数字的指数的对数

Question

我在做概率计算。我有很多非常小的数字，我想从1中减去所有的，然后精确地计算。我能准确地计算出这些小数字的对数。到目前为止，我的策略是这样的（使用numpy）：

给定一个小数字的对数数组x，计算：

y = numpy.logaddexp.reduce(x)

现在我想计算一些类似1-exp(y)或者更好的log(1-exp(y))，但是我不确定如何在不损失所有精度的情况下这样做。在

事实上，即使是logaddexp函数也遇到了精度问题。向量x中的值可以是-2到-800，甚至是更负的。上面的向量y基本上是在1e-16左右有一整段数字，这就是数据类型的eps。例如，精确计算的数据可以如下所示：

^{pr2}$

然后我试着计算指数的对数和：

In [359]: np.logaddexp.accumulate(x)
Out[359]: 
array([ -5.21946762e+00,  -3.66710221e+00,  -2.68983273e+00,
        -2.00815067e+00,  -1.51126604e+00,  -1.14067818e+00,
        -8.60829425e-01,  -6.48188808e-01,  -4.86276416e-01,
        -3.63085873e-01,  -2.69624488e-01,  -1.99028599e-01,
        -1.45996863e-01,  -1.06408884e-01,  -7.70565672e-02,
        -5.54467248e-02,  -3.96506186e-02,  -2.81859503e-02,
        -1.99225261e-02,  -1.40061296e-02,  -9.79701394e-03,
        -6.82045164e-03,  -4.72733966e-03,  -3.26317960e-03,
        -2.24396350e-03,  -1.53767347e-03,  -1.05026994e-03,
        -7.15209142e-04,  -4.85690052e-04,  -3.28980607e-04,
        -2.22305294e-04,  -1.49890553e-04,  -1.00858788e-04,
        -6.77380054e-05,  -4.54139175e-05,  -3.03974537e-05,
        -2.03154477e-05,  -1.35581905e-05,  -9.03659252e-06,
        -6.01552344e-06,  -3.99984336e-06,  -2.65671945e-06,
        -1.76283376e-06,  -1.16860435e-06,  -7.73997496e-07,
        -5.12213574e-07,  -3.38706792e-07,  -2.23809375e-07,
        -1.47785898e-07,  -9.75226648e-08,  -6.43149957e-08,
        -4.23904687e-08,  -2.79246430e-08,  -1.83858489e-08,
        -1.20995365e-08,  -7.95892319e-09,  -5.23300609e-09,
        -3.43929670e-09,  -2.25953475e-09,  -1.48391255e-09,
        -9.74194956e-10,  -6.39351406e-10,  -4.19466218e-10,
        -2.75121795e-10,  -1.80397409e-10,  -1.18254918e-10,
        -7.74993004e-11,  -5.07775611e-11,  -3.32619009e-11,
        -2.17835737e-11,  -1.42634249e-11,  -9.33764336e-12,
        -6.11190167e-12,  -3.99989955e-12,  -2.61737204e-12,
        -1.71253165e-12,  -1.12043465e-12,  -7.33052079e-13,
        -4.79645919e-13,  -3.13905885e-13,  -2.05519681e-13,
        -1.34650094e-13,  -8.83173582e-14,  -5.80300378e-14,
        -3.82338678e-14,  -2.52963381e-14,  -1.68421145e-14,
        -1.13181549e-14,  -7.70918073e-15,  -5.35155125e-15,
        -3.81152630e-15,  -2.80565548e-15,  -2.14872312e-15,
        -1.71971577e-15,  -1.43957518e-15,  -1.25665732e-15,
        -1.13722927e-15,  -1.05925916e-15,  -1.00835857e-15,
        -9.75131524e-16,  -9.53442707e-16,  -9.39286186e-16,
        -9.30046550e-16,  -9.24016349e-16,  -9.20080954e-16,
        -9.17512772e-16,  -9.15836886e-16,  -9.14743318e-16,
        -9.14029759e-16,  -9.13564174e-16,  -9.13260398e-16,
        -9.13062204e-16,  -9.12932898e-16,  -9.12848539e-16,
        -9.12793505e-16,  -9.12757603e-16,  -9.12734183e-16,
        -9.12718905e-16,  -9.12708939e-16,  -9.12702438e-16,
        -9.12698198e-16,  -9.12695432e-16,  -9.12693627e-16,
        -9.12692451e-16,  -9.12691683e-16,  -9.12691183e-16,
        -9.12690856e-16,  -9.12690643e-16,  -9.12690504e-16,
        -9.12690414e-16,  -9.12690355e-16,  -9.12690316e-16,
        -9.12690291e-16,  -9.12690275e-16,  -9.12690264e-16,
        -9.12690257e-16,  -9.12690252e-16,  -9.12690249e-16,
        -9.12690248e-16,  -9.12690246e-16,  -9.12690245e-16,
        -9.12690245e-16,  -9.12690245e-16,  -9.12690244e-16,
        -9.12690244e-16,  -9.12690244e-16,  -9.12690244e-16,
        -9.12690244e-16,  -9.12690244e-16,  -9.12690244e-16,
        -9.12690244e-16,  -9.12690244e-16,  -9.12690244e-16,
        -9.12690244e-16,  -9.12690244e-16,  -9.12690244e-16,
        -9.12690244e-16,  -9.12690244e-16,  -9.12690244e-16,
        -9.12690244e-16,  -9.12690244e-16,  -9.12690244e-16,
        -9.12690244e-16,  -9.12690244e-16,  -9.12690244e-16,
        -9.12690244e-16,  -9.12690244e-16,  -9.12690244e-16,
        -9.12690244e-16,  -9.12690244e-16,  -9.12690244e-16,
        -9.12690244e-16,  -9.12690244e-16,  -9.12690244e-16,
        -9.12690244e-16,  -9.12690244e-16,  -9.12690244e-16,
        -9.12690244e-16,  -9.12690244e-16,  -9.12690244e-16,
        -9.12690244e-16,  -9.12690244e-16,  -9.12690244e-16,
        -9.12690244e-16,  -9.12690244e-16,  -9.12690244e-16,
        -9.12690244e-16,  -9.12690244e-16,  -9.12690244e-16,
        -9.12690244e-16,  -9.12690244e-16,  -9.12690244e-16,
        -9.12690244e-16,  -9.12690244e-16,  -9.12690244e-16,
        -9.12690244e-16,  -9.12690244e-16,  -9.12690244e-16,
        -9.12690244e-16,  -9.12690244e-16,  -9.12690244e-16])

最终导致：

In [360]: np.logaddexp.reduce(x)
Out[360]: -9.1269024387687033e-16

所以我的精确性已经被抹去了。有什么办法解决这个问题吗？在

Answer 1

我对python了解不多，我的大部分工作都是用Java完成的，但在我看来，您最好先对所有值同时执行log sum exp技巧，而不是使用numpy.logaddexp.cumulate。在

在google中快速搜索可以在python库中找到一个候选库：scipy.misc.logsumexp。在

在任何情况下，编程都不难：

logsumexp(probs) := max(probs) + log(sum[i](exp(probes[i]-max(probs))))

比如说：

^{pr2}$

返回的单个值，比如p，只是probs中所有概率总和的对数。请注意，如果输入概率的最大值和较小值之间存在较大的差异，那么小的概率将被忽略，但是如果它们的贡献与大数相比非常小，那么在大多数应用程序中，它们的贡献应该是很小的。在

你可以使用更复杂的策略来计算这些小值，以防有大量的小数字，比如probe=0.5+1E-7+1E-7……数百万次，加起来就是0.1。你可以做的是把单个和分成几个数个，在合并之前先将大致相同比例的值相加。或者可以使用some-in-between pivot值而不是max，但必须确保最大值不太大，因为exp（probs[i]-pivot）在这种情况下会导致溢出。在

完成此操作后，您仍然需要进行calculatelog（1-exp（p））

为此，我发现这篇文章描述了一种使用逻辑函数尽可能避免精度损失的方法，这种方法可以在大多数语言的数学通用库中找到。在

Maechler M, Accurately Computing log(1 − exp(− |a|)) Assessed by the Rmpfr package, 2012

关键是根据输入值a使用两种可能的方法之一。在

定义：

log1mexp(a) := log(1-exp(a)) ### function that we seek to implement.
log1p(a) := log(1+a) # function provided by common math libraries.
exp1m(a) := exp(a) - 1 # function provided by common math libraries.

使用log1p实现log1mexp有一个明显的方法：

log1mexp(a) := log1p(-exp(a))

使用exp1m，您可以：

log1mexp(a) := log(-expm1(a))

然后，当a<；log（.5）时应使用log1p方法，当a>；=log（.5）时，应使用expm1方法。在

log1mexp(a) := (a < log(.5)) ? log1p(-exp(a)) : log(-expm1(a)).

更多信息请参考外部链接。在

Answer 2

我建议用它们的Taylor series替换exp（）和log（），相应地在0和1附近。这样，你就不会因为使用大数而失去精度（我刚才把1称为一个大数字：^）。使用Lagrange remainder公式或仅使用成员的表达式，并保留一些保留，以确定差异何时会超出您的精度。在

更新：

Python2.7的math.expm1（exp(x)-1）和math.log1p（log(1+x)）如果平台的C库的精度（通常是两倍）就足够了。（如果你不能用特殊的软件来计算FPU（如果你不能用特殊的软件来计算的话）。在

Answer 3

在Python2.7中，我们为这个用例添加了math.expm1()：

>>> from math import exp, expm1
>>> exp(1e-5) - 1  # gives result accurate to 11 places
1.0000050000069649e-05
>>> expm1(1e-5)    # result accurate to full precision
1.0000050000166668e-05

此外，在不损失精度的情况下，求和步骤还有math.fsum()：

^{pr2}$

最后，如果没有其他帮助，可以使用支持超高精度算术的decimal module：

>>> from decimal import *
>>> getcontext().prec = 200
>>> (1 - 1 / Decimal(7000000)).ln()
Decimal('-1.4285715306122546161332083855139723669559469615692284955124609122046580004888309867906750714454869716398919778588515625689415322789136206397998627088895481989036005482451668027002380442299229191323673E-7')