Python中的对数计算

Question

我正在寻找类似的计算：

其中f(i) 是一个函数，它在[-1,1] 中为{1,2,...,5000} 中的任何i 返回一个实数。

显然，总和的结果在 [-1,1] 的某个位置，但是当我似乎无法使用 Python 使用直接编码计算它时，0.5<sup>5000</sup> 变为 0 和 comb(5000,2000) 变为inf，这导致计算的总和变成NaN。

需要的解决方案是两边都使用log。

这是使用身份a × b = 2<sup>log(a) + log(b)</sup>，如果我可以计算log(a) 和log(b) 我可以计算总和，即使a 很大并且b 几乎是0。

所以我想我要问的是是否有一种简单的计算方法

log2(scipy.misc.comb(5000,2000))

所以我可以简单地计算我的总和

sum([2**(log2comb(5000,i)-5000) * f(i) for i in range(1,5000) ])

---

@abarnert 的解决方案在为 5000 数字工作时通过提高计算梳子的精度来解决问题。这适用于本示例，但无法扩展，因为如果我们使用 1e7 而不是 5000，则所需的内存将显着增加。

目前，我正在使用一种丑陋的解决方法，但可以保持低内存消耗：

log2(comb(5000,2000)) = sum([log2 (x) for x in 1:5000])-sum([log2 (x) for x in 1:2000])-sum([log2 (x) for x in 1:3000])

有没有办法在可读的表达式中做到这一点？

Answer 1

总和

是f 对binomial distribution 与n = 5000 和p = 0.5 的期望。

你可以用scipy.stats.binom.expect计算这个：

import scipy.stats as stats

def f(i):
    return i
n, p = 5000, 0.5
print(stats.binom.expect(f, (n, p), lb=0, ub=n))
# 2499.99999997

还请注意，当n 趋于无穷大，p 固定时，二项分布approaches the normal distribution 具有均值np 和方差np*(1-p)。因此，对于较大的n，您可以改为计算：

import math
print(stats.norm.expect(f, loc=n*p, scale=math.sqrt((n*p*(1-p))), lb=0, ub=n))
# 2500.0

Answer 2

编辑：@unutbu 已经回答了 real 的问题，但我会留在这里，以防 log2comb(n, k) 对任何人有用。

---

comb(n, k) 是 n！ / ((n-k)!k!) 和 n!可以使用Gamma functiongamma(n+1) 计算。 Scipy 提供了函数scipy.special.gamma。 Scipy还提供了gammaln，即Gamma函数的对数（即自然对数）。

所以log(comb(n, k)) 可以计算为gammaln(n+1) - gammaln(n-k+1) - gammaln(k+1)

例如log(comb(100, 8))（执行from scipy.special import gammaln后）：

In [26]: log(comb(100, 8))
Out[26]: 25.949484949043022

In [27]: gammaln(101) - gammaln(93) - gammaln(9)
Out[27]: 25.949484949042962

和日志(comb(5000, 2000)):

In [28]: log(comb(5000, 2000))  # Overflow!
Out[28]: inf

In [29]: gammaln(5001) - gammaln(3001) - gammaln(2001)
Out[29]: 3360.5943053174142

（当然，要获得以 2 为底的对数，只需除以 log(2)。）

为方便起见，您可以定义：

from math import log
from scipy.special import gammaln

def log2comb(n, k):
    return (gammaln(n+1) - gammaln(n-k+1) - gammaln(k+1)) / log(2) 

Answer 3

默认情况下，comb 给你一个float64，它溢出并给你inf。

但是，如果您通过 exact=True，它会为您提供一个 Python 可变大小的 int，它不会溢出（除非您变得如此可笑，以至于内存不足）。

而且，虽然您不能在 int 上使用 np.log2，但您可以使用 Python 的 math.log2。

所以：

math.log2(scipy.misc.comb(5000, 2000, exact=True))

---

作为替代方案，您的相对 n 选择 k 被定义为n!k / k!，对吗？你可以把它减少到∏(i=1...k)((n+1-i)/i)，这很容易计算。

或者，如果你想避免溢出，你可以交替使用* (n-i) 和/ (k-i)。

当然，您也可以简化为加减对数。我认为在 Python 中循环并计算 4000 个对数会比在 C 中循环并计算 4000 个乘法要慢，但我们总是可以对其进行向量化，然后可能会更快。让我们编写并测试：

In [1327]: n, k = 5000, 2000
In [1328]: %timeit math.log2(scipy.misc.comb(5000, 2000, exact=True))
100 loops, best of 3: 1.6 ms per loop
In [1329]: %timeit np.log2(np.arange(n-k+1, n+1)).sum() - np.log2(np.arange(1, k+1)).sum()
10000 loops, best of 3: 91.1 µs per loop

当然，如果你更关心记忆而不是时间……嗯，这显然会使情况变得更糟。我们一次有 2000 个 8 字节浮点数，而不是一个 608 字节整数。如果你增加到 100000、20000，你会得到 20000 个 8 字节浮点数，而不是一个 9K 整数。在 1000000、200000 时，它是 200000 个 8 字节浮点数与一个 720K 整数。

我不确定为什么这两种方式对您来说都是个问题。特别是考虑到您使用的是 listcomp 而不是genexpr，因此创建了一个不必要的 5000、100000 或 1000000 Python 浮点数列表——24MB 不是问题，但 720K 是？但如果是这样，我们显然可以迭代地做同样的事情，但会牺牲一些速度：

r = sum(math.log2(n-i) - math.log2(k-i) for i in range(n-k))

这并没有太比scipy 解决方案慢很多，而且它从不使用超过少量常量字节（少量 Python 浮点数）。 （除非您使用的是 Python 2，在这种情况下……只需使用 xrange 而不是 range，它就会恢复为常量。）

---

作为旁注，为什么是您使用列表推导而不是具有矢量化操作的 NumPy 数组（为了速度和紧凑性）或生成器表达式而不是列表推导（完全不使用内存，不以任何速度为代价）？