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