CODEWITHSUNDEEP
pythonscipyquad
mortysporty
mortysporty

Reputation: 2889

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

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

Function to be integrated <code>tpx_wx</code> over integration range

Upvotes: 2

Views: 1019

Answers (1)

user6655984
user6655984

Reputation:

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:

linux

On Windows:

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.

Upvotes: 1

Related Questions