scipy.integrate.quad fails (sometimes) when function to be integrated is also an integral

Question

I am using sqipy.integrate.quad to calculate a double integral. Basically I'm trying to calculate the integral over exp[-mu_wx_par] where mu_wx_par is also a integral.

My code mostly works. However, for some values it fails, i.e. it returns incorrect values.

import numpy as np
from scipy import integrate


def mu_wx_par(x, year, par):
    """ First function to be integrated """
    m = max(par['alfa'], 0) + par['beta'] * 10 ** (par['gamma'] * x)
    w = np.minimum(par['frem_a'] + par['frem_b'] * x + par['frem_c'] * x**2, 0)
    return m * (1 + w/100)**(year - par['innf_aar'])


def tpx_wx(x, t, par, year):
    """ Second function to be integrated (which contains an integral itself)"""
    result, _ = integrate.quad(lambda s: mu_wx_par(x + s, year + s, par), 0, t)
    return np.exp(-result)


def est_lifetime(x, year, par):
    """ Integral of second function. """
    result, _ = integrate.quad(lambda s: tpx_wx(x, s, par, year), 0, 125 - x)

    return result


# Test variables
par = {'alfa': 0.00019244401470947973,
       'beta': 2.420260552210541e-06,
       'gamma': 0.0525500987420195,
       'frem_a': 0.3244611019518985,
       'frem_b': -0.12382978382606026,
       'frem_c': 0.0011901237463116591,
       'innf_aar': 2018
       }


year = 2018
estimate_42 = est_lifetime(42, year, par)
estimate_43 = est_lifetime(43, year, par)
rough_estimate_42 = sum([tpx_wx(42, t, par, year) for t in range(0, 100)])

print(estimate_42)
print(estimate_43)
print(rough_estimate_42)
3.1184634065887544
46.25925442287578
47.86323490659588

The value of estimate_42 is not correct. It should be about the same value as rough_estimate_42. Note however that estimate_43 looks fine. What is going on here?

I am using scipy v1.1.0 and numpy v1.15.1 and Windows.

It has been suggested that the function is close to zero almost everywhere, as in this post scipy integrate.quad return an incorrect value. That is not the case as a simple plot of tpx_wx for x=42 from a=0 to b=125-42 clearly shows

from matplotlib import pyplot as plt

plt.plot(range(125-42), [tpx_wx(42, t, par, year) for t in range(0, 125-42)])
plt.show()

WIth SciPy 1.1.0, I get 47.36317494751913 for `estimate_42`, — , Oct 25 '18 at 13:59
Hmm. Im currently not at the office so I can't check my version og scipy. But I found that when I tried to do the integral for multiple x's (say for x between 20 and 80) some would fail whilst others would be ok. — mortysporty, Oct 25 '18 at 14:25
It is not the same issue. If that was the case, you would not have gotten the correct answer. There is nothing weird about the function at all. The fact that you get a different answer from me is suspiscious in itself. Running my code on a different computer also gave the incorrect answer. :| I use Windows btw. dont know if it matters. — mortysporty, Oct 26 '18 at 22:24
I also got 3.118... when on Windows. Strange. Anyway, adding `points=[(125-x)/2]` (in effect, bisecting the range of integration) fixed the issue. — , Oct 26 '18 at 23:58
Bisecting does not fix it. It merely makes quad avoid, in a totally random and unpredictable way, whathever numerical quirk or bug that caused the error. But anyway I appreciate your efforts. Thanks. — mortysporty, Oct 27 '18 at 06:41

score 1 · Accepted Answer · edited Jun 20 '20 at 09:12

This appears to be a known issue with the way that some Fortran code behind quad is compiled for Windows: calling it recursively can lead to failure in some cases. See also Big integration error with integrate.nquad.

Barring recompiling SciPy with better flags, it seems that one should avoid nested integration with quad while on Windows. One workaround is to use romberg method for one of the integration steps. Replacing quad in est_lifetime with

integrate.romberg(lambda s: tpx_wx(x, s, par, year), 0, 125 - x, divmax=20)

results in 47.3631754795 for estimate_42, consistent with what quad returns on Linux.

One way to visualize the process of integration is to declare a global list eval_points and insert eval_points.append(t) into tpx_wx. With the same version of SciPy (0.19.1 in this test), the results of plt.plot(eval_points, '.') look different.

On Linux:

On Windows:

Iterated bisection of the neighborhood of a tricky point around 60 is terminated prematurely on Windows, and it seems that the result thrown up is some partial integral over a subinterval.

Wow! You really went the extra mile here. Thank you. Using a different routine on one of the integrals is probably the way to go. Once again thanks! — mortysporty, Oct 28 '18 at 18:24

scipy.integrate.quad fails (sometimes) when function to be integrated is also an integral

1 Answers1