Python中的日志计算

Question

我想要计算类似的东西:

哪个f(i)函数返回[-1,1]任何iin中的实数{1,2,...,5000}.

显然,总和的结果在某处[-1,1],但是当我似乎无法使用直接编码在Python中计算它时,变为和变为,这导致计算的总和变成.0.550000comb(5000,2000)infNaN

所需的解决方案是使用双面登录.

这是使用身份,如果我可以计算,我可以计算总和,即使很大,几乎.a × b = 2log(a) + log(b)log(a)log(b)ab0

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

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

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

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

---

@ abarnert的解决方案,同时为5000图工作,通过提高计算梳子的精度来解决问题.这适用于此示例,但不会扩展,因为如果不是5000而我们需要1e7,所需的内存将显着增加.

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

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相对于一个二项式分布与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固定的情况下,二项分布接近具有均值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 已经回答了真正的问题，但我将把它留在这里，以防log2comb(n, k)对任何人有用。

comb(n, k)是n！/ ((nk)!k!) 和 n! 可以使用Gamma 函数 计算gamma(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

和 log(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)