Я использую 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 возвращает неверное значение. Это не так, поскольку простой график 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()
Хм. В настоящее время меня нет в офисе, поэтому я не могу проверить свою версию. Но я обнаружил, что, когда я пытался сделать интеграл для нескольких x (скажем, для x от 20 до 80), некоторые из них терпели неудачу, в то время как другие были в порядке.
... и это казалось довольно случайным, что не удастся.
Это не та же проблема. Если бы это было так, вы бы не получили правильный ответ. В этой функции нет ничего странного. Тот факт, что вы получили от меня другой ответ, сам по себе подозрительный. Выполнение моего кода на другом компьютере также дало неверный ответ. : | Кстати, я использую Windows. не знаю, имеет ли это значение.
Я также получил 3.118 ... когда на Windows. Странный. В любом случае, добавление points=[(125-x)/2]
(фактически, разделение диапазона интеграции пополам) устранило проблему.
Биссектриса не исправляет. Это просто заставляет quad совершенно случайным и непредсказуемым образом избегать числовой причуды или ошибки, вызвавшей ошибку. Но в любом случае я ценю ваши старания. Спасибо.
Это похоже на известная проблема из-за того, что некоторый код Fortran, лежащий в основе quad
, компилируется для Windows: его рекурсивный вызов в некоторых случаях может привести к сбою. См. Также Большая ошибка интеграции с Integrate.nquad.
Запрещая перекомпиляцию SciPy с лучшими флагами, кажется, что следует избегать вложенной интеграции с quad
в Windows. Один из обходных путей - использовать метод 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:
Итеративное деление пополам окрестности сложной точки около 60 преждевременно прерывается в Windows, и кажется, что полученный результат представляет собой некоторый частичный интеграл по подинтервалу.
Ух ты! Вы действительно прошли здесь лишнюю милю. Спасибо. Возможно, лучше использовать другую процедуру для одного из интегралов. Еще раз спасибо!
С SciPy 1.1.0 я получаю 47.36317494751913 для
estimate_42
,