scipy.integrate.quad 当要集成的函数也是一个整数时(有时)失败
scipy.integrate.quad fails (sometimes) when function to be integrated is also an integral
我正在使用 sqipy.integrate.quad
计算二重积分。基本上我试图计算 exp[-mu_wx_par] 的积分,其中 mu_wx_par 也是一个积分。
我的代码大部分都有效。但是,对于某些值它会失败,即 returns 不正确的值。
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
estimate_42
的值不正确。它应该与 rough_estimate_42
的值大致相同。但是请注意 estimate_43
看起来不错。这是怎么回事?
我正在使用 scipy v1.1.0 和 numpy v1.15.1 以及 Windows。
有人建议该函数几乎在所有地方都接近于零,如 post 中一样。事实并非如此,因为从 a=0
到 b=125-42
的 tpx_wx
for x=42
的简单图清楚地显示了
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()
这似乎是 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。
除非使用更好的标志重新编译 SciPy,似乎应该避免在 Windows 上与 quad
嵌套集成。一种解决方法是对其中一个集成步骤使用 romberg
方法。将 est_lifetime
中的 quad
替换为
integrate.romberg(lambda s: tpx_wx(x, s, par, year), 0, 125 - x, divmax=20)
对于 estimate_42
的 47.3631754795
结果与 quad
returns 在 Linux 上的结果一致。
可视化集成过程的一种方法是声明一个全局列表 eval_points
并将 eval_points.append(t)
插入 tpx_wx
。使用相同版本的SciPy(本次测试为0.19.1),plt.plot(eval_points, '.')
的结果看起来不同。
在 Linux:
在 Windows:
围绕 60 的棘手点邻域的迭代二分在 Windows 上过早终止,似乎抛出的结果是子区间上的部分积分。
我正在使用 sqipy.integrate.quad
计算二重积分。基本上我试图计算 exp[-mu_wx_par] 的积分,其中 mu_wx_par 也是一个积分。
我的代码大部分都有效。但是,对于某些值它会失败,即 returns 不正确的值。
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
estimate_42
的值不正确。它应该与 rough_estimate_42
的值大致相同。但是请注意 estimate_43
看起来不错。这是怎么回事?
我正在使用 scipy v1.1.0 和 numpy v1.15.1 以及 Windows。
有人建议该函数几乎在所有地方都接近于零,如 post a=0
到 b=125-42
的 tpx_wx
for x=42
的简单图清楚地显示了
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()
这似乎是 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。
除非使用更好的标志重新编译 SciPy,似乎应该避免在 Windows 上与 quad
嵌套集成。一种解决方法是对其中一个集成步骤使用 romberg
方法。将 est_lifetime
中的 quad
替换为
integrate.romberg(lambda s: tpx_wx(x, s, par, year), 0, 125 - x, divmax=20)
对于 estimate_42
的 47.3631754795
结果与 quad
returns 在 Linux 上的结果一致。
可视化集成过程的一种方法是声明一个全局列表 eval_points
并将 eval_points.append(t)
插入 tpx_wx
。使用相同版本的SciPy(本次测试为0.19.1),plt.plot(eval_points, '.')
的结果看起来不同。
在 Linux:
在 Windows:
围绕 60 的棘手点邻域的迭代二分在 Windows 上过早终止,似乎抛出的结果是子区间上的部分积分。