我正在使用sqipy.integrate.quad计算双整数。基本上，我试图在exp [-mu_wx_par]上计算积分，其中mu_wx_par也是积分。
我的代码大部分有效。但是，对于某些值，它会失败，即返回错误的值。
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。
有人建议该函数几乎在所有地方都接近于零，如这篇文章scipy integrate.quad return an incorrect value。情况并非如此，因为
tpx_wx从x=42到a=0的b=125-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，是为Windows编译quad后面的某些Fortran代码的方式：在某些情况下递归调用它可能导致失败。另请参见Big integration error with integrate.nquad。
除非使用更好的标志来重新编译SciPy，否则似乎应该避免在Windows上与quad嵌套集成。一种解决方法是对集成步骤之一使用romberg方法。将quad中的est_lifetime替换为
integrate.romberg(lambda s: tpx_wx(x, s, par, year), 0, 125 - x, divmax=20)
导致
47.3631754795的estimate_42，与quad在Linux上返回的一致。可视化集成过程的一种方法是声明全局列表
eval_points并将eval_points.append(t)插入tpx_wx。使用相同版本的SciPy（本测试中为0.19.1），plt.plot(eval_points, '.')的结果看起来有所不同。在Linux上：
在Windows上：
Windows上过早地终止了一个棘手点附近的迭代二等分，大约在60终止，看来在一个子区间上抛出的结果是部分积分。
关于python - scipy.integrate.quad有时(有时)要集成的功能也是不可或缺的，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/52983963/