python - 尝试使用 scipy quad 集成函数并使用 brentq 求解时出现数学范围错误

Question

我希望 scipy 为我计算以下积分：

def F(x, theta, mu, sigma, r):
    #equation 1
    def integrand(u):
        return u**(r/mu - 1) * exp(sqrt(2*mu / sigma**2) * (x-theta)*u - u**2/2)
    return scipy.integrate.quad(integrand, 1e-5, np.inf)[0]

def Prime(f, x, theta, mu, sigma, r, h=1e-5):
    # given f, estimates f'(x) using the difference quotient formula 
    return (f(x+h, theta, mu, sigma, r) - f(x, theta, mu, sigma, r)) / h

def b_star(theta, mu, sigma, r, c):
    #equation 2
    def func(b):
        return F(b, theta, mu, sigma, r) - (b-c)*Prime(F, b, theta, mu, sigma, r)
    
    return scipy.optimize.brentq(func, -2, 3)

对应于以下等式 ，

其中我为 b 求解方程 2。这似乎在大多数情况下都有效，但是在某些情况下会发生以下错误：

In: b_star(0.573085,3.213728,0.137655,0,0.001)
Out:
---------------------------------------------------------------------------
OverflowError                             Traceback (most recent call last)
<ipython-input-164-1aeb31cb3219> in <module>
----> 1 b_star(0.573085,3.213728,0.137655,0,0.001)

<ipython-input-163-14ce7e008b20> in b_star(theta, mu, sigma, r, c)
     38         return F(b, theta, mu, sigma, r) - (b-c)*Prime(F, b, theta, mu, sigma, r)
     39 
---> 40     return so.brentq(func, -2, 3)
     41 
     42 def V(x, theta, mu, sigma, r, c):

/Library/Frameworks/Python.framework/Versions/3.8/lib/python3.8/site-packages/scipy/optimize/zeros.py in brentq(f, a, b, args, xtol, rtol, maxiter, full_output, disp)
    778     if rtol < _rtol:
    779         raise ValueError("rtol too small (%g < %g)" % (rtol, _rtol))
--> 780     r = _zeros._brentq(f, a, b, xtol, rtol, maxiter, args, full_output, disp)
    781     return results_c(full_output, r)
    782 

<ipython-input-163-14ce7e008b20> in func(b)
     36 
     37     def func(b):
---> 38         return F(b, theta, mu, sigma, r) - (b-c)*Prime(F, b, theta, mu, sigma, r)
     39 
     40     return so.brentq(func, -2, 3)

<ipython-input-163-14ce7e008b20> in F(x, theta, mu, sigma, r)
     18     def integrand(u):
     19         return u**(r/mu - 1) * exp(sqrt(2*mu / sigma**2) * (x-theta)*u - u**2/2)
---> 20     return si.quad(integrand, 1e-5, np.inf)[0]
     21 
     22 def G(x, theta, mu, sigma, r):

/Library/Frameworks/Python.framework/Versions/3.8/lib/python3.8/site-packages/scipy/integrate/quadpack.py in quad(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst)
    339 
    340     if weight is None:
--> 341         retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit,
    342                        points)
    343     else:

/Library/Frameworks/Python.framework/Versions/3.8/lib/python3.8/site-packages/scipy/integrate/quadpack.py in _quad(func, a, b, args, full_output, epsabs, epsrel, limit, points)
    453             return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
    454         else:
--> 455             return _quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit)
    456     else:
    457         if infbounds != 0:

<ipython-input-163-14ce7e008b20> in integrand(u)
     17     # equation 3.3
     18     def integrand(u):
---> 19         return u**(r/mu - 1) * exp(sqrt(2*mu / sigma**2) * (x-theta)*u - u**2/2)
     20     return si.quad(integrand, 1e-5, np.inf)[0]
     21 

OverflowError: math range error

我尝试了几种解决方案，例如使用 Decimal 包或 gmpy 包来允许更大的浮点数，但在这些情况下输出变为 0.0。我还尝试将输入转换为 np.float128 数据类型，这也没有解决问题。

有解决办法吗？

标签： python