Fixed the issue by adding a fallback to bissection

[Cedriq1astaken]

close #2061

Issue

The issue described by jherekhealy occurs in the quantile_newton method when dealing with extreme value near 0 or 1 (more often near 1). It fails to converge to a value due to floating point error and the default tolerance being too small resulting in either an infinite loop as it oscillates between values, or converging very slowly.

Fix

I modified quantile_newton, cquantile_newton, invlogcdf_newton, invlogccdf_newton by adding:

• a new parameter maxiter (which by default = 10_000), which defines a maximum iteration before it falls back
• A cycle detector that detects whenever a 2-cycles occurs

When either cases occur, quantile_newton falls back to quantile_bisect

Test

I wrote a Monte Carlo simulation to test extreme values near 1, [0.9999; 0.999999999999]

Here are the results for 1,000,000 samples at tol = 1e-12:

quantile_newton (no fallback)
  converged              859955   86.00%
  oscillation            140045   14.00%

quantile_newton (with bisection fallback)
  converged              1000000  100.00%

The test shows that the fallback correctly detects the 2-cycles and converges successfully every time instead of infinitely looping.

AI Assistance Disclosure

The debugging strategy, proposed fix, and interpretation of the numerical issue were developed by the author. An AI coding tool was used to help draft and refine parts of the implementation, tests, and PR text. All AI-assisted changes were reviewed, edited, and validated by the author before submission.

[codecov-commenter]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 86.75%. Comparing base (03ddda7) to head (7820dc3).
⚠️ Report is 2 commits behind head on master.

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #2066      +/-   ##
==========================================
+ Coverage   86.63%   86.75%   +0.11%     
==========================================
  Files         149      149              
  Lines        8921     8887      -34     
==========================================
- Hits         7729     7710      -19     
+ Misses       1192     1177      -15     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[andreasnoack]

The solution here is too crude. The more fundamental problem is that the stopping rule doesn't take into account the granularity of the floating point numbers around the solution. Try to ask the AI to consider that and about how this is traditionally handled in quantile functions. I also think we should have a max_inter, but it should be something like 50 or 100, and if reached we should at least warn and possibly error. Ideally, you'd never hit it.

[Cedriq1astaken]

Hi, thanks for the comment,
I've done what you suggested and the AI suggested in response to scale stopping rule with the magnitude of x,
and add a granularity floor using cbrt(eps)^2 to avoid chasing precision we can't achieve.

Also dropped max_iter to 50 and added a fallback to bisection with warning when things don't converge cleanly.
Is that alright?

[devmotion]

I think we should use Roots.jl. This would also easily allow us to use better bracketing root solvers such as ITP instead of vanilla bisection.

[Cedriq1astaken]

Okay that makes sense., I switched the fallback to use Roots.jl with ITP instead.

[devmotion]

My point was to generally make use of Roots. Being able to use ITP instead of vanilla bisection is just one benefit. Generally, the algorithms and numerical corner cases should be much better battle tested in Roots.

[Cedriq1astaken]

Just to make sure I understand: are you suggesting the we replace the internal implementation of quantile_newton to just call Roots.jl instead of using the custom loop?"

[devmotion]

Yes, that's my suggestion. We would tune only the options if the default ones are not appropriate.

[Cedriq1astaken]

Okay, I've now refactored the core logic of quantile_newton and the others to use Roots.jl and I upgraded the quantile_bisect to use ITP too.

[devmotion]

Before going through the code, can you add all currently failing quantile-related examples, such as #2061, #1611, #1571, #1999, #1807, #1869 (maybe there are more), as tests? Then it would be clearer whether the fixes are sufficient and we would avoid silent regressions in these cases.

[Cedriq1astaken]

I added regression tests for the six quantile-related issues you listed. I found a couple of other quantile-related issues like #1415, but they don't seem related to the convergence issues this PR focuses on. Should I fix and include those too?

[devmotion]

Completed by #2070