求幂过程中的溢出错误答案

【问题标题】：OverFlowError during exponentiation求幂过程中的溢出错误
【发布时间】：2017-11-01 13:24:39
【问题描述】：

我无法计算这个。该代码适用于所有 N 值，包括 N = 57，但对于大于或等于 58 的 N 会引发溢出错误（34，“结果太大”）。有没有办法解决这个问题？谢谢。

import numpy as np
import scipy.integrate
import scipy.optimize
import warnings
warnings.filterwarnings('ignore')

R = 6
N = 58
Nb = 4
S_norm = 0.3
Ao = 1/8.02e-5


y = (N-Nb/R)/Ao

def likelihood(psi):

    def func_likelihood(A):
        sum = 0
        for k in range(0, N+1):
            r = (np.math.factorial(N-k+Nb)/(np.math.factorial(k)*np.math.factorial(N-k))*(A*psi*(1+R))**k)
            sum = r+sum
        return sum* (np.exp(-0.5*(np.log(A/Ao)/S_norm)**2)*np.exp(-A*psi)/A)

    return (scipy.integrate.quad(func_likelihood, 0, np.inf, full_output=1,epsabs=1e-300, epsrel=1e-300))[0]

psi = y-y/10

MLE = scipy.optimize.fmin(lambda psi: -likelihood(psi), y, disp=0,full_output=1)[0][0]

normal_factor = 1/likelihood(MLE)

print(normal_factor* likelihood(psi))

【问题讨论】：

整数除法可能会有所帮助。例如，(np.math.factorial(N-k+Nb) // (np.math.factorial(k)*np.math.factorial(N-k)) 或 (np.math.factorial(N-k+Nb) // (np.math.factorial(k) // np.math.factorial(N-k))。顺便说一句，您不应该使用 sum 作为变量名，因为这会影响内置的 sum 函数。
在 scipy 或 numpy 中还应该有一个优化的二项式函数。
@PM2Ring，整数除法不会给出正确的结果。感谢您关于使用“总和”作为变量的建议。
哦，好吧，这是一个尝试。 ;) 您是否尝试重新排列这些术语，例如(factorial(N-k+Nb) // factorial(N-k)) / factorial(k)？另一种选择是在 range(N-k+Nb, N-k, -1) 的循环中执行大乘法，然后找到该结果的 GCD 和 factorial(k)，这样您就可以在尝试浮点除法之前用整数除法消除公因数。
请注意，np.math.factorial 的拼写有点误导，因为它与 NumPy 无关：它只是常规的标准库 math.factorial。 (np.math is math -> True)

标签： python math integer-overflow exponentiation

【解决方案1】：

计算项的商并使用该商在每个步骤中更新项通常更有效且溢出问题更少。

压缩索引k的术语是

r[k] = ( (N-k+Nb)! )/( k!*(N-k)! ) * (A*psi*(1+R))**k

所以最后一项的商是

r[k+1] / r[k] = ( (N-k) )/( (N-k+Nb)*(k+1) ) * (A*psi*(1+R))

和

r[0] = ( (N+Nb)! )/( N! ) = (N+1)*...*(N+Nb)

然后允许将求和函数重新表述为

 def func_likelihood(A):
      sum = 0
      r = 1
      for k in range(Nb): r *= N+k+1
      sum = r
      for k in range(0, N+1):
          sum += r;
          r *= (A*psi*(1+R)) * (N-k)/((N-k+Nb)*(k+1))
      return sum * (np.exp(-0.5*(np.log(A/Ao)/S_norm)**2)*np.exp(-A*psi)/A)

【讨论】：

【解决方案2】：

您未能提供错误消息，但我打赌它是在从 for k 循环内的阶乘中计算 r。

有几种方法可以解决这个问题。一种是切换运算顺序，这样你就消除了公因数，而不是从58开始！（超出默认整数限制）。这涉及更多编码，或者可能调用组合阶乘函数，例如将执行 C(n, k) 的函数——@LutzL 提到了二项式函数。当然，您并没有相当使用二项式，但您可以使用它，然后根据需要调整分子以适应 -k+Nb 因子。

另一种方法是切换到大整数（请参阅 numpy 和 scipy 的文档，以便您解决所有包中的问题）。

【讨论】：

np.math.factorial 已经使用了大整数，希望这些阶乘可以安全地划分。但我同意使用正确的C(n, k) 更好。
也许——但不是处理 58！是一个可疑的截止日期。
是的，但那是因为在大整数上使用/ 会产生浮点结果。