求幂期间的 OverFlowError

Question

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

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))

Answer 1

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

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

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

Answer 2

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

索引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)

求幂期间的 OverFlowError

OverFlowError during exponentiation

python

math

integer-overflow

exponentiation