math range error when trying to integrate a function using scipy quad and solving using brentq

Question

I'd like SciPy to solve the following problem for me:

import numpy as np
from numpy import exp, sqrt
import scipy.integrate
import scipy.optimize

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_star(0.573085, 3.213728, 0.5, 0, 0.001)  # 0.7714361802263615

This corresponds with the following equations:

,

in which I solve the last equation for b.

This seems to work most of the time; however, in some circumstances the following error occurs:

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

I tried several solutions like using the Decimal package or the gmpy package to allow for bigger floating point numbers, but the output becomes 0.0 in those cases.

I also tried converting the inputs to a np.float128 datatype, but this didn't solve the problem either.

Is there a fix for this?

Answer 1

As @BoarGules notes, the argument to exp is too large and the result overflows.

Running your code verbatim except for the replacement:

# ~scipy.optimize.brentq(func, -2, 3)~ F(3) overflows
scipy.optimize.brentq(func, -2, 2.5)

reveals that there is a root near 0.5636522268679528. This is not bigger than theta, but you can check for yourself that it is a root of the provided equation.

If you can't tighten the bracket you provide to brentq, you can use mpmath to avoid overflow, replacing calls to scipy.integrate.quad and scipy.optimize.brentq with mpmath's quad and findroot.

import scipy
import numpy as np
from mpmath import exp, sqrt, log, mp
mp.dps = 30
# not all `mpf` conversions are strictly necessary, but
# sometimes it's better to be safe.
mpf = mp.mpf

def F(x, theta, mu, sigma, r):
    #equation 1
    def integrand(u):
        return u**(r/mu - mp.one) * exp(sqrt(2*mu / sigma**2) * (x-theta)*u - u**2/2)

    return mp.quad(integrand, (mpf(1e-5), mp.inf))

def Prime(f, x, theta, mu, sigma, r, h=1e-5):
    h = mpf(h)
    # 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):
    theta, mu, sigma, r, c = mpf(theta), mpf(mu), mpf(sigma), mpf(r), mpf(c)
    #equation 2
    def func(b):
        return F(b, theta, mu, sigma, r) - (b-c)*Prime(F, b, theta, mu, sigma, r)
    
    return mp.findroot(func, x0=[mpf(-2), mpf(3)], solver='bisect')

b_star(0.573085, 3.213728, 0.137655, 0, 0.001)

You can also replace your implementation of Prime with use of mpmath's diff for a more accurate derivative:

def Prime(f, x, theta, mu, sigma, r):
    return mp.diff(lambda x: f(x, theta, mu, sigma, r), x)

There are ways of performing this calculation in log-space with SciPy, but the mpmath solution is more straightforward.