Fast Newton-Raphson inverse_cdf for Chi and InverseGamma#404
Conversation
Chi and InverseGamma were the only continuous distributions still using the generic bisection default for inverse_cdf. Give them a custom solver in the same brent-like + Newton-Raphson vein as Gamma (statrs-dev#382): a shared internal::newton_raphson_quantile that brackets the quantile to a factor of two and refines it with safeguarded Newton steps, falling back to bisection whenever a step is non-finite or leaves the bracket. Two refinements over a plain port keep it accurate and fast across the whole range, including the tails statrs-dev#390 cared about: - convergence is tested on the relative step, not prec::convergence's absolute 1e-9 (meaningless for a 1e-12 quantile); and the Newton step is checked for convergence before the bracket is tightened, so a converged step that rounds onto an endpoint is not rejected into a spurious bisection. - the upper half inverts sf rather than cdf, which saturates to one and loses the resolution to place a deep upper-tail quantile. Matches scipy/mpmath (dps=60) to <= 3.4e-15 relative error over p in [1e-12, 1 - 1e-12] for both distributions, and converges in a handful of Newton steps (medians ~5, vs ~48 fixed bisection steps), roughly 60% fewer cdf/pdf evaluations across the grid.
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #404 +/- ##
==========================================
+ Coverage 94.77% 94.82% +0.04%
==========================================
Files 61 61
Lines 13402 13539 +137
==========================================
+ Hits 12702 12838 +136
- Misses 700 701 +1 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
The safeguarded step rejected a Newton iterate only when it was non-finite or outside the bracket, never when it was simply not making progress. That is enough for a quadratically converging tail but not for a linear one: on the inverse gamma's exp(-rate/x) lower tail each step advances rate/x by exactly one however far the root is, so the search crept toward the quantile and ran out of iterations short of it — InverseGamma(1, 1).inverse_cdf(1e-200) came back 4.3% high, and InverseGamma(1, 400).inverse_cdf(1e-300) 9.0% high, both of which the generic bisection this replaces got right. Restore Numerical Recipes' second rtsafe test: take the Newton step only while it also shrinks at least as fast as a bisection would. The bracket then halves at least every other iteration, so the budget always suffices. Nine deep-lower-tail cases across four parameter sets now agree with mpmath (dps=60) to under 2e-16 relative. Also test the p == 0, p == 1 and out-of-range arms of both new inverse_cdf implementations, which had no coverage.
|
Chasing the uncovered lines turned up a real defect rather than a coverage gap: the safeguarded step rejected a Newton iterate only when it was non-finite or left the bracket, never when it simply was not making progress, so on the inverse gamma's exp(-rate/x) lower tail (each step advances rate/x by one however far the root is) the search ran out of iterations short of the root — InverseGamma(1, 1).inverse_cdf(1e-200) came back 4.3% high and InverseGamma(1, 400).inverse_cdf(1e-300) 9.0% high, both of which the bisection this replaces got right. Restoring rtsafe's second test, take the Newton step only while it also shrinks at least as fast as a bisection, fixes those; nine deep-tail cases across four parameter sets now match mpmath at dps=60 to under 2e-16, and with tests added for the p == 0, p == 1 and out-of-range arms patch coverage is 100%. |
1 participant
Follows #390 / #382, per @YeungOnion's suggestion to extend the brent-like + Newton-Raphson
inverse_cdffrom Gamma (#382) to the distributions that still fall back to the generic default.ChiandInverseGammawere the last two continuous distributions using the generic bisection default forinverse_cdf(the two named in #390). This gives each a custom solver in the same vein as Gamma.Approach
A shared
internal::newton_raphson_quantile(p, cdf, sf, pdf)brackets the quantile to a factor of two, then refines it with safeguarded Newton steps (Numerical Recipes'rtsafe):[low, high]is only ever tightened, and a Newton step that is non-finite or leaves the bracket falls back to bisection.ChiandInverseGammacall it with their own cdf/sf/pdf.Two things were needed to keep it accurate and fast across the whole range, not just the middle:
prec::convergencefolds in an absolute1e-9tolerance, which is meaningless when the quantile itself is ~1e-12, so the solver tests the relative Newton step instead. It also checks the raw Newton step for convergence before tightening the bracket — otherwise a converged step that rounds onto the endpoint just moved toxis rejected and the search bisects the whole span back (this showed up as ~30 wasted iterations near the median).sfin the upper half. Ascdfsaturates to one it can no longer resolve a deep upper-tail quantile (InverseGamma's reaches ~6e23atp = 1 − 1e-12); invertingsfthere keeps it well conditioned, as in Converge the default continuous inverse_cdf instead of a fixed 16 iterations #390.Validation
Checked against
scipy.stats, cross-checked with mpmath at 60 digits, overp ∈ [1e-12, 1 − 1e-12]forChi(k), k ∈ {1,2,3,5,10} andInverseGamma(a,b), (a,b) ∈ {(1,1),(0.5,0.5),(2,1),(3,2),(5,3)}:3.4e-15; round-tripcdf(inverse_cdf(p)) ≈ pto ~3e-15.The
test_inverse_cdf_reference/test_inverse_cdf_round_tripcases added in #390 (scipy/mpmath values pinned in both tails) now exercise the new solver and still pass; a newinternal::test_newton_raphson_quantilecovers the shared routine directly against the Exponential closed form.cargo test,cargo fmt --checkandcargo clippyare clean.Gamma keeps its own #382 solver here; happy to fold it into this helper too as a follow-up if you'd prefer a single routine.