Skip to content

Converge the default continuous inverse_cdf instead of a fixed 16 iterations - #390

Merged
YeungOnion merged 1 commit into
statrs-dev:mainfrom
gaoflow:fix-continuous-default-inverse-cdf-convergence
Jul 19, 2026
Merged

Converge the default continuous inverse_cdf instead of a fixed 16 iterations#390
YeungOnion merged 1 commit into
statrs-dev:mainfrom
gaoflow:fix-continuous-default-inverse-cdf-convergence

Conversation

@gaoflow

@gaoflow gaoflow commented Jul 18, 2026

Copy link
Copy Markdown
Contributor

The ContinuousCDF::inverse_cdf default runs a fixed 16-iteration bisection with no
convergence test, so it only ever narrows the initial bracket by 2^16 regardless of how
wide the bracket is or how small the quantile is. Every continuous distribution without a
closed-form quantile inherits this — currently Chi and InverseGamma — and the result
sits far from the crate's own DEFAULT_RELATIVE_ACC (1e-14):

  • Chi::new(1).inverse_cdf(1e-12) returns 3.05e-5 instead of 1.25e-12; in fact every
    p <= ~1e-6 returns the same 3.05e-5 bracket floor. Even the median is only accurate
    to ~1e-5.
  • InverseGamma shows the same ~1e-5 baseline error, growing into the tails.

This is the continuous mirror of the discrete side, which already converges via
internal::integral_bisection_search.

The fix replaces the fixed count with a bisection that runs until the bracket agrees to
DEFAULT_RELATIVE_ACC (capped at 100 iterations and by float resolution). Two smaller
changes make the tails accurate:

  • Bracket from the distribution's own min()/max() when finite, only doubling out from
    ±2 when a bound is infinite. Chi/InverseGamma are supported on [0, ∞), so the old
    [-2, 2] seed wasted the lower half.
  • In the upper half, invert sf rather than cdf. As cdf saturates to one it can no
    longer place the quantile (many x map to a single f64 cdf value); sf stays well
    conditioned. Distributions that don't override sf fall back to 1 - cdf, so this is a
    no-op for them.

Across p ∈ [1e-12, 1-1e-12] the worst relative error for Chi and InverseGamma drops
from ~1e7 (and ~1e-5 at the median) to <= 1.3e-14. Reference quantiles in the new tests
come from scipy, cross-checked against mpmath at 60 digits; a round-trip test covers a
wider grid.

let q = Chi::new(1).unwrap().inverse_cdf(1e-12); // was 3.05e-5, want 1.25e-12

@codecov

codecov Bot commented Jul 18, 2026

Copy link
Copy Markdown

Codecov Report

❌ Patch coverage is 99.05660% with 1 line in your changes missing coverage. Please review.
✅ Project coverage is 94.94%. Comparing base (7f6c9ee) to head (49983c2).
⚠️ Report is 10 commits behind head on main.

Files with missing lines Patch % Lines
src/distribution/mod.rs 97.95% 1 Missing ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##             main     #390      +/-   ##
==========================================
+ Coverage   94.73%   94.94%   +0.21%     
==========================================
  Files          59       59              
  Lines       13052    13146      +94     
==========================================
+ Hits        12365    12482     +117     
+ Misses        687      664      -23     

☔ 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.

The ContinuousCDF::inverse_cdf default ran a fixed 16-iteration bisection
with no convergence check, narrowing the bracket by only 2^16 regardless of
its width. Distributions without a closed-form quantile (Chi, InverseGamma)
inherited errors far above the crate's DEFAULT_RELATIVE_ACC:
Chi::new(1).inverse_cdf(1e-12) returned 3.05e-5 instead of 1.25e-12, and even
the median was only accurate to ~1e-5.

Bisect until the bracket meets DEFAULT_RELATIVE_ACC, seed it from the
distribution's finite domain bounds, and invert sf in the upper half where cdf
saturates to one. Worst-case relative error over p in [1e-12, 1-1e-12] drops
from ~1e7 to <=1.3e-14 for both affected distributions.
@gaoflow
gaoflow force-pushed the fix-continuous-default-inverse-cdf-convergence branch from 3c07a7c to 49983c2 Compare July 18, 2026 15:04
@YeungOnion

Copy link
Copy Markdown
Contributor

Thanks!

would you be willing follow this with something similar to the brent-like approach that uses NR steps that you introduced for Gamma in #382?

@YeungOnion
YeungOnion merged commit 5738f71 into statrs-dev:main Jul 19, 2026
11 checks passed
YeungOnion pushed a commit that referenced this pull request Jul 27, 2026
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.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants