我正在尝试计算 (3e28 选择 2e28)/2^(3e28)。我试过 scipy.misc.comb 来计算 3e28 选择 2e28 但它给了我 inf。当我计算 2^(3e28) 时，它引发了 OverflowError: (34, 'Result too large')。如何计算或估计 (3e28 选择 2e28)/2^(3e28)？
【问题讨论】：
标签： python numbers scipy calculator
以下使用我的回答here中的log2comb：
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)
log2p = log2comb(3e28, 2e28) - 3e28
print "log2p =", log2p
打印出来的
log2p = -2.45112497837e+27
因此，您的数字的以 2 为底的 对数 约为 -2.45e27。如果您尝试计算 2**log2p，您会得到 0。也就是说，该数字小于标准 64 位浮点数可表示的最小正数。
【讨论】：
您可以使用大n 的二项式的正态近似来计算此比率。当n 很大时，k 必须相对接近n/2 才能使(n choose k) / 2^n 不可忽略。
这里有一些代码可以计算这个：
def n_choose_k_over_2_pow_n(n, k):
# compute the mean and standard deviation of the normal
# approximation
mu = n / 2.
sigma = np.sqrt(n) * 1/4.
# now transform to a standard normal variable
z = (k - mu) / sigma
return 1/np.sqrt(2*np.pi) * np.exp(-1/2. * z**2)
这样：
>>> n_choose_k_over_2_pow_n(3e28, 2e28)
0.0
>>> n_choose_k_over_2_pow_n(3e28, 1.5e28)
0.3989422804014327
如您所见，计算下溢。一个解决方案是计算答案的log，我们可以用这段代码来做：
def log_n_choose_k_over_2_pow_n(n, k):
# compute the mean and standard deviation of the normal
# approximation
mu = n / 2.
sigma = np.sqrt(n) * 1/4.
# now transform to a standard normal variable
z = (k - mu) / sigma
# return the log of the answer
return -1./2 * (np.log(2 * np.pi) + z**2)
另一个快速检查：
>>> log_n_choose_k_over_2_pow_n(3e28, 2e28)
-6.6666666666666638e+27
>>> log_n_choose_k_over_2_pow_n(3e28, 1.5e28)
-0.91893853320467267
如果我们对这些取幂，我们将得到之前的答案。
我们可以通过求助于统计结果来做到这一点。二项分布由下式给出：
P(K = k) = (n choose k) p^k * p^(n-k)
对于大的n，这很好地近似于正态分布，均值n*p，方差n*p*(1-p)。
将p 设置为1/2。然后我们有：
P(K = k) = (n choose k) (1/2)^k * (1/2)^(n-k)
= (n choose k) (1/2)^n
= (n choose k) / (2^n)
这正是你的比率的形式。因此，在转换为均值 n/2 和方差 n/4 的标准正态变量后，我们可以通过对标准正态分布 pdf 的简单评估来计算您的比率。
【讨论】：
(3e28)^-1。
e^(-1.699e+27)，这绝对是微不足道的。这比 (3e28)^-1（大约 e^-65.571）小了 10^20 个数量级。 （这些数字有时真的让我头疼……）
(3e28 choose (3e28/2))/(2^(3e28))。根据您的推理，它“基本上为零”。但事实上并非如此。 （我在计算它时遇到了麻烦，因为我在日志空间中得到了灾难性的取消......这通常表明结果在日志空间中接近于零，这意味着它们不可忽略）
有一些 Python 库可以让您进行任意精度的算术运算。例如 SymPy 中使用的 mpmath。
不过，您必须重写代码才能使用库函数。
http://docs.sympy.org/latest/modules/mpmath/basics.html?highlight=precision
编辑：我刚刚注意到您正在处理的数字的大小 - 对于我的建议来说太大了。
【讨论】：
使用斯特林近似（在 1e10+ 范围内非常准确），结合对数：
(3e28 choose 2e28) / 2^(3e28) = 3e28! / [(3e28 - 2e28)! * 2e28!] / 2^(3e28)
= e^ [log (3e28!) - log((3e28-2e28)!) - log(2e28!) - 3e28 * log(2)]
并从那里应用斯特林的近似值：
log n! ~= log(sqrt(2*pi*n)) + n*log(n) - n
你会得到答案的。
下面是这个近似值准确度的示例：
>>> import math
>>> math.log(math.factorial(100))
363.73937555556347
>>> math.log((2*math.pi*100)**.5) + 100*math.log(100) - 100
363.7385422250079
对于 100!，它在日志空间中下降了不到 0.01%。
【讨论】：