改善与涉及集成

我尝试使用下面的代码（不相关的部分去掉）来解决积分方程方程scipy.optimize.fsolve精度：改善与涉及集成

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使得从tau到无穷大的pdf的积分等于flag。请注意，x0是脚本中其他位置定义的根（浮点）估计值。还请注意，flag是一个非常小的数字，大小为1e-9。

在我的应用程序中，fsolve只能成功找到约50％的时间。它通常只会返回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)得到一个断言错误。我不太确定如何解决这个问题。任何帮助将非常感谢！

来源

2017-03-22 oirectine

如果我正确地读取了您的代码，您设置它的方式会浪费9位数的精度;难怪scipy正在挣扎。你不能从_t_到+ infinity执行积分吗？ –

@PaulPanzer对不起，但9位数字从哪里来？ – oirectine

如果我正确理解你，你最终希望匹配_flag_，这是1e-9的顺序。但是你在1 - _flag_上做了实际的根发现，这意味着如果你用相对精度为1e-12（这是相当苛刻的并且不总是可能的）然后从1 - _S_变回到_S_，然后_S_将会有1e-12阶的绝对误差，所以相对误差为1e-3阶。如果你可以从_tau_到+ infinity集成，那么你可以直接与_flag_比较，这个问题就会消失。 –

举个例子，我想1/x为一体1和alpha之间检索目标整体2.0。这

import quadpy 
from scipy.optimize import fsolve 

def f(alpha): 
    beta, _ = quadpy.line_segment.adaptive_integrate(
     lambda x: 1.0/x, [1, alpha], 1.0e-10 
     ) 
    return beta 

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

正确返回7.3890561。

的失败quadpy断言

assert all(lengths > minimum_interval_length)

你得到手段，自适应集成打了极限：要么放松你的宽容一点，或减少minimum_interval_length（见here）。

来源

2017-04-10 15:17:25

谢谢！我从此得到这个与其他集成方法一起工作，但想尝试并获得quadpy工作 - 我已经对上面的帖子进行了编辑。 – oirectine

@oirectine改编我的答案。 –

再次感谢 - 我已经试过放松宽容和减少minimum_interval_length无济于事。我对自适应正交不太了解;做某些功能可能不适合它？我也尝试在具有已知值的'fsolve'之外运行'adaptive_integrate'函数，并且得到了一些无意义的结果。任何帮助深表感谢！ – oirectine

改善与涉及集成

回答

相关问题