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

Question

I'm trying to solve an integral equation using the following code (irrelevant parts removed):

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]

Which takes a probability density function and finds a result tau such that the integral of pdf from tau to infinity is equal to flag. Note that x0 is a (float) estimate of the root defined elsewhere in the script. Also note that flag is an extremely small number, on the order of 1e-9.

In my application fsolve only successfully finds a root about 50% of the time. It often just returns x0, significantly biasing my results. There is no closed form for the integral of pdf, so I am forced to integrate numerically and feel that this might be introducing some inaccuracy?

EDIT:

This has since been solved using a method other than that described below, but I'd like to get quadpy to work and see if the results improve at all. The specific code I'm trying to get to work is as follows:

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)

When I run this I get an assertion error from assert all(lengths > minimum_interval_length). I'm not quite sure of how to remedy this; any help would be very much appreciated!

Answer 1

As an example, I tried 1 / x for the integration between 1 and alpha to retrieve the target integral 2.0. This

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)

correctly returns 7.38905611.

The failing quadpy assertion

assert all(lengths > minimum_interval_length)

you're getting means that the adaptive integration hit its limit: Either relax your tolerance a bit, or decrease the minimum_interval_length (see here).