BUG: Inaccuracy in `logcdf` with `quadrature` method

[JacobHass8]

I was working on writing more tests and the hypothesis library found this example where the logcdf method wasn't accurate with the quadrature method. Here is the script

import numpy as np
from normal import Normal

if __name__ == '__main__':
    arg = -0.55686482
    methods = ['formula', 'log/exp', 'quadrature']
    X = Normal(mu=0.904775196502295, sigma=1.3102019638947042)
    reference = np.exp(X.logcdf(arg))

    print("-------------------")
    for method in methods:
        result = np.exp(X.logcdf(arg, method=method))
        print(f"Method: {method}", "\n"
              "Absolute difference: ", np.abs(reference - result), "\n"
              "-------------------")

and the output is

-------------------
Method: formula 
Absolute difference:  0.0 
-------------------
Method: log/exp 
Absolute difference:  1.1102230246251565e-16 
-------------------
Method: quadrature 
Absolute difference:  1.3001154173353235e-05 
-------------------

The relative tolerance is currently set to 1e-7 so this doesn't pass. It's not as if this is in the tail so I don't understand why the accuracy isn't that great.

[mdhaber]

Hypothesis is very good at finding these sorts of counterexamples. It doesn't indicate an actionable bug in the integrator or implementation¹; the integrator is just not accurate to machine precision with default settings for arbitrary arguments. There's not much we can do about it besides adjust the tolerance. This is one of the reasons we need to move away from CI testing with hypothesis - occasional reminders of this imperfection puts pressure to improve, but it's not something that should hold up CI for everyone.

Footnotes

1. I have studied this particular problem in the past and found there are little regions where the tanh-sinh rule just doesn't integrate the normal PDF well with a low number of points, but it so happens that doubling the number of points once doesn't change the answer much. This makes it think the answer is accurate, so it terminates. If we tighten the tolerance or otherwise force the number of points to be higher, the problem goes away, but it would be unfortunate to make everything slower across the board for these unusual cases. ↩

[mdhaber]

That said, this does motivate some of the improvements to the generic tests noted in #1 (comment) : this is exactly a case of "the particular numbers used are very unlucky, and there's nothing we can reasonably change about the implementation to avoid the failure", so we need to find a way to nudge the test just a little bit to avoid the uninformative failure without otherwise weakening (e.g. skipping, adding huge tolerances) the test.