From c3fe60ff64d796c86af725e81dafcb539f89097a Mon Sep 17 00:00:00 2001 From: Vincent Gao Date: Mon, 20 Jul 2026 09:50:15 +0200 Subject: [PATCH 1/2] Add safeguarded-Newton inverse_cdf for Chi and InverseGamma 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 (#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 #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. --- src/distribution/chi.rs | 26 ++++++++ src/distribution/internal.rs | 107 ++++++++++++++++++++++++++++++ src/distribution/inverse_gamma.rs | 27 ++++++++ 3 files changed, 160 insertions(+) diff --git a/src/distribution/chi.rs b/src/distribution/chi.rs index c9482bc6..6ec26bad 100644 --- a/src/distribution/chi.rs +++ b/src/distribution/chi.rs @@ -148,6 +148,32 @@ impl ContinuousCDF for Chi { gamma::gamma_ur(self.freedom() as f64 / 2.0, x * x / 2.0) } } + + /// Calculates the inverse cumulative distribution function for the chi + /// distribution at `p`, i.e. the `p`-quantile. + /// + /// # Panics + /// + /// If `p` is not in `[0, 1]`. + fn inverse_cdf(&self, p: f64) -> f64 { + if !(0.0..=1.0).contains(&p) { + panic!("p must be in [0, 1]") + } + if p == 0.0 { + return self.min(); + } + if p == 1.0 { + return self.max(); + } + // The chi cdf has no closed-form inverse; solve it with the shared + // safeguarded Newton search instead of the generic bisection default. + super::internal::newton_raphson_quantile( + p, + |x| self.cdf(x), + |x| self.sf(x), + |x| self.pdf(x), + ) + } } impl Min for Chi { diff --git a/src/distribution/internal.rs b/src/distribution/internal.rs index 9b146c37..3cb6c793 100644 --- a/src/distribution/internal.rs +++ b/src/distribution/internal.rs @@ -37,6 +37,81 @@ pub fn integral_bisection_search( } } +/// Quantile `F^{-1}(p)` for a continuous distribution supported on `(0, ∞)` +/// whose cdf has no closed-form inverse, found with a safeguarded Newton–Raphson +/// step (Numerical Recipes' `rtsafe`); `cdf`, `sf` and `pdf` are the +/// distribution's own functions and `p` must lie strictly inside `(0, 1)`. +/// +/// The bracket `[low, high]` is kept as an invariant and only ever tightened, so +/// a Newton step that is non-finite or would leave the bracket falls back to +/// bisection and the iterate can never escape the support — the guard that keeps +/// a quantile sitting against a boundary from diverging to NaN (cf. Gamma in +/// #382). In the upper half the survival function is inverted rather than the +/// cdf: as `cdf` saturates to one it can no longer resolve a deep upper-tail +/// quantile, whereas `sf` stays well conditioned there. +pub fn newton_raphson_quantile( + p: f64, + cdf: impl Fn(f64) -> f64, + sf: impl Fn(f64) -> f64, + pdf: impl Fn(f64) -> f64, +) -> f64 { + // Bracket the quantile within a factor of two, `cdf(low) <= p <= cdf(high)`, + // by walking a unit interval out toward it. Moving *both* ends keeps the + // bracket tight when the quantile is far from 1 (deep in either tail), so the + // Newton phase starts close and converges in a handful of steps rather than + // bisecting the whole span. + let mut low = 1.0; + let mut high = 2.0; + while cdf(low) > p { + high = low; + low /= 2.0; + } + while cdf(high) < p { + low = high; + high *= 2.0; + } + + // Solve `sf(x) = 1 - p` in the upper half and `cdf(x) = p` otherwise; either + // way the residual is increasing in `x` with derivative `pdf(x)`. + let upper = p > 0.5; + let target = if upper { 1.0 - p } else { p }; + + // A *relative* accuracy target: quantiles span many orders of magnitude, so + // an absolute tolerance (as in `prec::convergence`) is meaningless deep in a + // tail where the quantile itself is far smaller than that tolerance. + let accuracy = crate::prec::DEFAULT_RELATIVE_ACC; + const MAX_ITERATIONS: usize = 100; + let mut x = (low + high) / 2.0; + for _ in 0..MAX_ITERATIONS { + let residual = if upper { + target - sf(x) + } else { + cdf(x) - target + }; + let newton = x - residual / pdf(x); + // A full Newton step below the relative tolerance means we have + // converged; accept it before the bracket bookkeeping below, which would + // otherwise reject a converged step that has rounded onto the very + // endpoint we are about to move to `x`. + if (newton - x).abs() <= accuracy * x.abs() { + return newton; + } + // Tighten the bracket by the sign of the (increasing) residual, then take + // the Newton step only while it stays strictly inside, else bisect. + if residual >= 0.0 { + high = x; + } else { + low = x; + } + x = if newton.is_finite() && newton > low && newton < high { + newton + } else { + (low + high) / 2.0 + }; + } + x +} + #[cfg(test)] macro_rules! testing_boiler { ($($arg_name:ident: $arg_ty:ty),+; $dist:ty; $dist_err:ty) => { @@ -460,6 +535,38 @@ mod test { } } + #[cfg(feature = "std")] + #[test] + fn test_newton_raphson_quantile() { + use super::newton_raphson_quantile; + // Exponential(λ) has the closed-form quantile -ln(1 - p)/λ, so it is an + // exact oracle for the shared solver across both tails. + for lambda in [0.5f64, 1.0, 4.0] { + let cdf = |x: f64| -(-lambda * x).exp_m1(); + let sf = |x: f64| (-lambda * x).exp(); + let pdf = |x: f64| lambda * (-lambda * x).exp(); + for p in [ + 1e-12, + 1e-8, + 1e-4, + 1e-2, + 0.25, + 0.5, + 0.75, + 0.99, + 1.0 - 1e-6, + 1.0 - 1e-12, + ] { + let got = newton_raphson_quantile(p, cdf, sf, pdf); + let want = -(-p).ln_1p() / lambda; // -ln(1 - p)/λ, accurate in both tails + assert!( + (got - want).abs() <= 1e-12 * want, + "Exp({lambda}).inverse_cdf({p}) = {got}, want {want}" + ); + } + } + } + pub mod boiler_tests { use crate::distribution::{Beta, BetaError}; use crate::statistics::*; diff --git a/src/distribution/inverse_gamma.rs b/src/distribution/inverse_gamma.rs index 7b231d0e..1f2a3faa 100644 --- a/src/distribution/inverse_gamma.rs +++ b/src/distribution/inverse_gamma.rs @@ -173,6 +173,33 @@ impl ContinuousCDF for InverseGamma { gamma::gamma_lr(self.shape, self.rate / x) } } + + /// Calculates the inverse cumulative distribution function for the inverse + /// gamma distribution at `p`, i.e. the `p`-quantile. + /// + /// # Panics + /// + /// If `p` is not in `[0, 1]`. + fn inverse_cdf(&self, p: f64) -> f64 { + if !(0.0..=1.0).contains(&p) { + panic!("p must be in [0, 1]") + } + if p == 0.0 { + return self.min(); + } + if p == 1.0 { + return self.max(); + } + // The inverse gamma cdf has no closed-form inverse; solve it with the + // shared safeguarded Newton search instead of the generic bisection + // default, which loses precision far into the (very heavy) upper tail. + super::internal::newton_raphson_quantile( + p, + |x| self.cdf(x), + |x| self.sf(x), + |x| self.pdf(x), + ) + } } impl Min for InverseGamma { From a0984e94a3cdafc7682d71bc2c3c6dd00b852b21 Mon Sep 17 00:00:00 2001 From: Vincent Gao Date: Sun, 26 Jul 2026 18:55:09 +0200 Subject: [PATCH 2/2] Bound the Newton iteration count and cover the quantile guards MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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. --- src/distribution/chi.rs | 19 ++++++++++++++ src/distribution/internal.rs | 36 +++++++++++++++++--------- src/distribution/inverse_gamma.rs | 42 +++++++++++++++++++++++++++++++ 3 files changed, 85 insertions(+), 12 deletions(-) diff --git a/src/distribution/chi.rs b/src/distribution/chi.rs index 6ec26bad..f4dd5abc 100644 --- a/src/distribution/chi.rs +++ b/src/distribution/chi.rs @@ -579,4 +579,23 @@ mod tests { } } } + + #[test] + fn test_inverse_cdf_p0_p1() { + let d = create_ok(3); + assert_eq!(d.inverse_cdf(0.0), d.min()); + assert_eq!(d.inverse_cdf(1.0), d.max()); + } + + #[test] + #[should_panic(expected = "p must be in [0, 1]")] + fn test_inverse_cdf_p_above_one() { + create_ok(3).inverse_cdf(1.0 + f64::EPSILON); + } + + #[test] + #[should_panic(expected = "p must be in [0, 1]")] + fn test_inverse_cdf_p_below_zero() { + create_ok(3).inverse_cdf(-1e-300); + } } diff --git a/src/distribution/internal.rs b/src/distribution/internal.rs index 3cb6c793..792b3b64 100644 --- a/src/distribution/internal.rs +++ b/src/distribution/internal.rs @@ -43,12 +43,14 @@ pub fn integral_bisection_search( /// distribution's own functions and `p` must lie strictly inside `(0, 1)`. /// /// The bracket `[low, high]` is kept as an invariant and only ever tightened, so -/// a Newton step that is non-finite or would leave the bracket falls back to -/// bisection and the iterate can never escape the support — the guard that keeps -/// a quantile sitting against a boundary from diverging to NaN (cf. Gamma in -/// #382). In the upper half the survival function is inverted rather than the -/// cdf: as `cdf` saturates to one it can no longer resolve a deep upper-tail -/// quantile, whereas `sf` stays well conditioned there. +/// a Newton step that is non-finite, would leave the bracket, or is not shrinking +/// at least as fast as a bisection falls back to bisection; the iterate can never +/// escape the support — the guard that keeps a quantile sitting against a +/// boundary from diverging to NaN (cf. Gamma in #382) — and the bracket is halved +/// often enough to reach the root within the iteration budget. In the upper half +/// the survival function is inverted rather than the cdf: as `cdf` saturates to +/// one it can no longer resolve a deep upper-tail quantile, whereas `sf` stays +/// well conditioned there. pub fn newton_raphson_quantile( p: f64, cdf: impl Fn(f64) -> f64, @@ -82,6 +84,7 @@ pub fn newton_raphson_quantile( let accuracy = crate::prec::DEFAULT_RELATIVE_ACC; const MAX_ITERATIONS: usize = 100; let mut x = (low + high) / 2.0; + let mut last_step = high - low; for _ in 0..MAX_ITERATIONS { let residual = if upper { target - sf(x) @@ -96,18 +99,27 @@ pub fn newton_raphson_quantile( if (newton - x).abs() <= accuracy * x.abs() { return newton; } - // Tighten the bracket by the sign of the (increasing) residual, then take - // the Newton step only while it stays strictly inside, else bisect. + // Tighten the bracket by the sign of the (increasing) residual. if residual >= 0.0 { high = x; } else { low = x; } - x = if newton.is_finite() && newton > low && newton < high { - newton + // Take the Newton step only while it stays strictly inside the bracket + // *and* shrinks at least as fast as a bisection would, else bisect. The + // second half of that test is what bounds the iteration count: a tail + // where the step is merely linear — the inverse gamma's `exp(-b/x)` + // advances `b/x` by one per step regardless of how far the root is — + // otherwise creeps toward the root and runs out of iterations short of + // it, e.g. `InverseGamma(1, 1).inverse_cdf(1e-200)` off by 4%. + let step = (newton - x).abs(); + if newton.is_finite() && newton > low && newton < high && 2.0 * step <= last_step { + x = newton; + last_step = step; } else { - (low + high) / 2.0 - }; + last_step = (high - low) / 2.0; + x = low + last_step; + } } x } diff --git a/src/distribution/inverse_gamma.rs b/src/distribution/inverse_gamma.rs index 1f2a3faa..b7ba4cf5 100644 --- a/src/distribution/inverse_gamma.rs +++ b/src/distribution/inverse_gamma.rs @@ -540,4 +540,46 @@ mod tests { } } } + + #[test] + fn test_inverse_cdf_deep_lower_tail() { + // A Newton step on `exp(-rate/x)` only advances `rate/x` by one however + // far the root is, so without the bisection-progress guard the search + // runs out of iterations short of it. References from mpmath (dps=60). + let cases: &[(f64, f64, f64, f64)] = &[ + (1.0, 1.0, 1e-100, 0.0043429448190325185), + (1.0, 1.0, 1e-200, 0.0021714724095162593), + (1.0, 1.0, 1e-300, 0.0014476482730108394), + (0.01, 1.0, 1e-200, 0.00222287682578071), + (0.01, 1.0, 1e-300, 0.0014711980613782147), + (100.0, 1.0, 1e-100, 0.0020694527796494867), + (100.0, 1.0, 1e-200, 0.0013193408747368574), + (1.0, 400.0, 1e-300, 0.5790593092043358), + (0.5, 1e6, 1e-200, 2188.7520668986626), + ]; + for &(a, b, p, expected) in cases { + let q = InverseGamma::new(a, b).unwrap().inverse_cdf(p); + let relerr = ((q - expected) / expected).abs(); + assert!(relerr <= 1e-14, "InverseGamma({a}, {b}).inverse_cdf({p}) = {q}, want {expected} (relerr {relerr:e})"); + } + } + + #[test] + fn test_inverse_cdf_p0_p1() { + let d = create_ok(1.0, 1.0); + assert_eq!(d.inverse_cdf(0.0), d.min()); + assert_eq!(d.inverse_cdf(1.0), d.max()); + } + + #[test] + #[should_panic(expected = "p must be in [0, 1]")] + fn test_inverse_cdf_p_above_one() { + create_ok(1.0, 1.0).inverse_cdf(1.0 + f64::EPSILON); + } + + #[test] + #[should_panic(expected = "p must be in [0, 1]")] + fn test_inverse_cdf_p_below_zero() { + create_ok(1.0, 1.0).inverse_cdf(-1e-300); + } }