使用涉及积分的方程式提高 scipy.optimize.fsolve 的准确性

Question

我正在尝试使用以下代码求解积分方程（已删除无关部分）：

def _pdf(self, a, b, c, t):
    pdf = some_pdf(a,b,c,t)
    return pdf

def _result(self, a, b, c, flag):
    return fsolve(lambda t: flag - 1 + quad(lambda tau: self._pdf(a, b, c, tau), 0, t)[0], x0)[0]

它采用概率密度函数并找到结果 tau，使得 pdf 从 tau 到无穷大的积分等于 flag。请注意，x0 是脚本中其他地方定义的根的（浮动）估计值。另请注意，flag 是一个非常小的数字，大约为 1e-9。

在我的应用程序中，fsolve 仅在大约 50% 的时间内成功找到根。它通常只是 returns x0，严重影响了我的结果。 pdf的积分没有封闭形式，所以我被迫进行数值积分，感觉这可能会引入一些不准确的地方？

编辑：

此后已使用下述方法以外的方法解决了此问题，但我想让 quadpy 工作并查看结果是否有所改善。我尝试开始工作的具体代码如下：

import quadpy
import numpy as np
from scipy.optimize import *
from scipy.special import gammaln, kv, gammaincinv, gamma
from scipy.integrate import quad, simps

l = 226.02453163
mu = 0.00212571582056
nu = 4.86569872444
flag = 2.5e-09
estimate = 3 * mu

def pdf(l, mu, nu, t):
    return np.exp(np.log(2) + (l + nu - 1 + 1) / 2 * np.log(l * nu / mu) + (l + nu - 1 - 1) / 2 * np.log(t) + np.log(
        kv(nu - l, 2 * np.sqrt(l * nu / mu * t))) - gammaln(l) - gammaln(nu))


def tail_cdf(l, mu, nu, tau):
    i, error = quadpy.line_segment.adaptive_integrate(
        lambda t: pdf(l, mu, nu, t), [tau, 10000], 1.0e-10
    )
    return i

result = fsolve(lambda tau: flag - tail_cdf(l, mu, nu, tau[0]), estimate)

当我运行这样做时，我从 assert all(lengths > minimum_interval_length) 得到一个断言错误。我不太确定如何解决这个问题；任何帮助将不胜感激！

Answer 1

作为示例，我尝试 1 / x 在 1 和 alpha 之间进行积分以检索目标积分 2.0。这个

import quadpy
from scipy.optimize import fsolve

def f(alpha):
    beta, _ = quadpy.quad(lambda x: 1.0/x, 1, alpha)
    return beta

target = 2.0
res = fsolve(lambda alpha: target - f(alpha), x0=2.0)
print(res)

正确 returns 7.38905611.

失败的 quadpy 断言

assert all(lengths > minimum_interval_length)

你得到的意思是自适应集成达到了极限：稍微放宽你的容忍度，或者降低 minimum_interval_length（参见 here）。

使用涉及积分的方程式提高 scipy.optimize.fsolve 的准确性

Improving accuracy in scipy.optimize.fsolve with equations involving integration

python

solver

scipy

numerical-methods

numerical-integration