改善與涉及集成

我嘗試使用下面的代碼（不相關的部分去掉）來解決積分方程方程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

改善與涉及集成

回答

相關問題