1

我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 integration.quad return an wrong value。情况并非如此，因为tpx_wxfor x=42from a=0to的简单情节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()

要在积分范围内积分的函数<code>tpx_wx</code>

4

1 回答 1

1

这似乎是为 Windows 编译某些 Fortran 代码的方式的一个已知问题quad：在某些情况下，递归调用它可能会导致失败。另请参阅Integrate.nquad 的大集成错误。

除非使用更好的标志重新编译 SciPy，否则似乎应该避免quad在 Windows 上与 while 进行嵌套集成。一种解决方法是romberg对集成步骤之一使用方法。替换quad为est_lifetime_

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

结果47.3631754795为 for ，与Linux 上的返回结果estimate_42一致。quad

可视化集成过程的一种方法是声明一个全局列表eval_points并插入eval_points.append(t)到tpx_wx. 使用相同版本的 SciPy（本次测试为 0.19.1），结果plt.plot(eval_points, '.')看起来不同。

在 Linux 上：

linux

在 Windows 上：

视窗

在 Windows 上，60 附近的一个棘手点的邻域的迭代二等分过早终止，似乎抛出的结果是某个子区间上的部分积分。

于 2018-10-27T21:37:12.677 回答