取幂期间的OverFlowError

问题描述：

我在计算这个时遇到了麻烦。该代码适用于所有的N值（包括N = 57），但是当N大于或等于58时会引发溢出错误（34，'Result too large'）。有没有办法解决这个问题？谢谢。取幂期间的OverFlowError

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（Nk））'或'（np.math.factorial（N-k + Nb）//（np.math.factorial（k）// np.math.factorial（Nk））'。顺便说一句，你不应该使用'sum'作为变量名称，因为这会影响内置的sum '函数 –

还应该有一个优化的二项式函数在scipy或numpy中可用 – LutzL

@ PM2Ring，整数除法并没有给出正确的结果感谢您使用'sum'作为变量的建议 – DPdl

答

它通常更高效，并且具有更少的溢出问题来计算术语的商并使用商来更新每个步骤中的术语。

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

这然后允许重新配制求和函数as

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)

这很好用！非常感谢LutzL！ :) – DPdl

答

你没有提供的错误信息，但我打赌它在for k循环中计算从阶乘r。

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

另一种方法是切换到大整数（请参阅文档以获得更多信息，以便您解决所有包中的问题）。

'np.math.factorial'已经使用大整数，希望这些因子可以安全地分开。但我同意使用适当的'C（n，k）'更好。 –

也许 - 但不处理58！是一个可疑的切断点。 – Prune

确实如此，但是这是因为在大整数上使用'/'会产生浮点结果。 –

取幂期间的OverFlowError

相关推荐