fix: use Loader's saddle-point form for the gamma/beta prefixes and pmfs#418
Open
agene0001 wants to merge 4 commits into
Open
fix: use Loader's saddle-point form for the gamma/beta prefixes and pmfs#418agene0001 wants to merge 4 commits into
agene0001 wants to merge 4 commits into
Conversation
`approx`'s `ulps_eq` short-circuits on `abs_diff_eq(epsilon)` before it
consults the ULPs distance. The crate's wrapper paired it with
`DEFAULT_EPS = 1e-9`, which made the ULPs bound unreachable, so every internal
`ulps_eq!(x, y)` was really a 1e-9 absolute comparison.
That matters because the macro is used to recognise *exact* parameter values, so
anything within 1e-9 of them silently took a degenerate branch:
Binomial(p = 1 - 2^-33, n = 100).pmf(99) 0 (true 1.16e-8)
Binomial(p = 1 - 2^-33, n = 100).ln_pmf(99) -inf
Geometric(p = 1 - 2^-33).max() 1 (true u64::MAX)
Geometric(p = 1 - 2^-33).skewness() inf (true 92681.9)
Beta(1 + 2^-33, 1 + 2^-33).pdf(0.0) 1 (true 0)
digamma(-1 + 2^-33) -inf (true -8.59e9)
Adds `DEFAULT_ULPS_EPS = f64::EPSILON` and uses it for `ulps_eq!` and
`assert_ulps_eq!`. The exact values still take the degenerate branches -
`Binomial(1.0, 5).pmf(5) == 1`, `Geometric(1.0).max() == 1`, `digamma` at the
poles is still -inf - all covered by existing tests.
Tests use `1 - 2^-33` rather than `1 - 1e-10` so that `1 - p` is exact and the
references are not limited by the representation of `p`.
Also fixes `Binomial::entropy`, which summed `-p * ln(p)` unguarded and so
returned `NaN` once any mass underflowed to zero (`0 * -inf`) - for `p = 0.5`
that is every `n` past about 1100. `Categorical::entropy` already filtered zero
terms; this brings `Binomial` in line.
`ln_gamma`/`gamma` evaluate Pugh's 10-term Lanczos approximation as the
partial-fraction sum `d_0 + sum_k d_k / (z + k - 1)`. The residues alternate in
sign, so that sum cancels badly - condition number ~3600 around z = 50, i.e.
~440 eps of relative error in the sum alone.
The same quantity as a single fraction `N(z) / D(z)`, derived from those residues
in exact rational arithmetic, has *all-positive* coefficients in both numerator
and denominator (`D(z) = z (z+1) ... (z+9)`, whose expanded coefficients are
unsigned Stirling numbers of the first kind and exact in f64). For z > 0 both
Horner evaluations then have condition number 1. This needs no new external
constants - they are derived from the residues already in the file - and agrees
with the partial-fraction form to 3.2e-17 over [0.5, 3000]. A reversed form in
1/z covers z >= 1e29, where z^10 would overflow.
Two further fixes in the same evaluation:
* `lanczos_power` compensates the base of `((p + g) / e)^p`. `powf` amplifies a
relative error in its base by the exponent, so the two roundings in
`(p + g) / e` were the dominant error (~190 ulp by x = 122). The residuals of
the addition and the division - the latter against a double-double `e` - are
recovered exactly and applied as a first-order correction, which is only
valid while it stays small, so it is gated on that.
* `gamma`/`ln_gamma` at the positive integers up to 171 come from the exact
factorial table, which also makes `Gamma(2) == 1` rather than
1.0000000000000002.
* `gamma` for x above ~169.7 halves the exponent and squares, because the power
alone overflows there while the full product stays representable to
x ~ 171.61. `gamma(171.6)` was `inf`; it is now 1.5858969096673e308.
`sin_pi`/`tan_pi` reduce by the period before calling `sin`/`tan`. `x - round(x)`
is exact, so the error stays relative to the fractional part instead of being
`ulp(PI) / dist_to_pole`, which near the poles cost several decimal digits.
Measured against mpmath at 45 dps:
before after
ln_gamma p99 45.9 ulp 8.2 ulp
gamma (positive) 285 median 19 median
gamma (negative) 109 median, 1173 max 3.4 median, 22.8 max
digamma (negative) 19.5 median, 1.5e10 1.6 median, 538 max
`gamma`'s remaining ~19-35 ulp is the approximation floor of the f64-rounded Pugh
coefficient set itself (~3.7e-15, confirmed by evaluating the formula in exact
arithmetic), so the evaluation is now compensated to well below the fit.
Four test expectations move as a consequence, each verified against mpmath: one
`Binomial::sf` and two `FisherSnedecor` pdf/ln_pdf literals were fitted to the
old output, and the `Dirichlet`/`MultivariateStudent` doctest values were 34 and
22 ulp from truth (the new outputs are within 2).
The Lentz recurrence in `checked_beta_reg` was capped at a fixed 140 iterations
and silently returned whatever it had reached. It is slowest at the centre of the
distribution, where the worst case over `x` grows like `5 * min(a, b).cbrt()`, so
past `min(a, b) ~ 1.5e4` the answer was simply wrong. Against the exact identity
`I_{1/2}(a, a) == 1/2`:
a = b = 1e5 0.49999969504 relative error 6.1e-7
a = b = 1e6 0.49121972700 1.8e-2
a = b = 1e7 0.21285001452 5.7e-1
`Binomial::new(0.5, 2e6).cdf(1e6)` returned 0.4916 against a true 0.50028.
The bound now scales as `8 * min(a, b).cbrt()`, measured to cover the worst case
over `x` with headroom, clamped to bound the work at roughly 10 ms. The loop still
exits on convergence, so ordinary calls are unaffected: `beta_reg(2.5, 2.5, 0.5)`
is 150 ns against 145 ns, and `beta_reg(1, 1, 0.3)` is unchanged.
Three robustness fixes alongside it, all found by probing the parameter space
rather than by sweeping accuracy:
* an underflowed prefix now short-circuits to the corresponding endpoint. The
result is `bt * h / a` with `h` of order one, so this is exact - and it avoids
forming `0.0 * h`, which is NaN whenever the recurrence overflowed
(`beta_reg(1e300, 1e-300, 0.5)` was NaN on both sides of this change before).
* `x == 0` and `x == 1` return 0 and 1 directly. They used to depend on the
symmetry test, which mapped `x == 0` to `1.0` once `a + b` overflowed.
* the symmetry threshold `(a + 1) / (a + b + 2)` is computed scaled when
`a + b` overflows, since otherwise it collapses to zero and sends every `x`
down the transformed branch.
* the result is kept in `[0, 1]`, and falls back to the concentrated-limit step
function if the truncated recurrence produced something non-finite. `I_x` is a
probability and callers such as `Binomial::cdf` are contractually so;
`a = b = 1e20` previously returned -2.56.
Regression tests use two exact identities that need no reference data -
`I_{1/2}(a, a) == 1/2` and `I_x(a, b) + I_{1-x}(b, a) == 1` - plus a 12x12x6
parameter grid asserting the result stays a finite probability. Both identity
tests fail on the old fixed bound.
Densities and incomplete-function prefixes in this family are written as
`exp(a ln x - x - ln_gamma(a))` and friends. Every term there grows like
`a ln a` while the sum stays `O(ln a)`, so at `a = 1e4` each term is ~1e5 and the
~1e-11 absolute error of the sum becomes the accuracy ceiling of everything
downstream. That is what limited `beta_reg` to ~6e-11 even where its recurrence
converged.
Adds the two building blocks from Loader, "Fast and Accurate Computation of
Binomial Probabilities" (2000):
* `stirling_delta(z)`, the Stirling-series remainder of `ln_gamma`, which is
`O(1/(12z))`. It is table-free: the recurrence
`d(z) = d(z+1) + (z + 1/2) ln(1 + 1/z) - 1` lifts any argument into the range
where a 6-term series is good to 1.4e-18, so there is no interpolation table
and no recursive dependency on `ln_gamma`. Accurate to ~1e-15 *absolute*,
which is what matters since every caller adds it to an exponent.
* `bd0(x, np) = x ln(x/np) - (x - np)`, which has a double root at `x == np`.
Evaluated by a series in `(x - np) / (x + np)` so the cancellation is done
analytically, giving ~1e-16 relative accuracy uniformly through the root.
Rewriting the prefixes in terms of these keeps every piece `O(1)`. Measured
against mpmath, and against `I_{1/2}(a, a) == 1/2` which is exact:
beta_reg, a = b = 1e4 6e-11 -> 8.7e-15
beta_reg, a = b = 1e6 3e-9 -> 2.6e-13
gamma_ur, a >= 16 5.3 median ulp -> 0.38, p99 970k -> 3071
Binomial::pmf 685k median ulp -> 2.45
Poisson::pmf 649 median ulp -> 3.16
Both prefixes fall back to the direct form for small parameters, where there is
no cancellation left to remove and `stirling_delta`'s recurrence would cost more
roundings than it saves - without that gate, `FisherSnedecor(1,1).cdf` got
*worse*, since a = b = 0.5 needs 16 recurrence steps.
`bd0` is sensitive enough to its mean argument that the rounding of `n * p` alone
matters: `d bd0 / d np = 1 - x / np`, so half an ulp of `np` at 6e5 moves the
result 2.5e-13, which was a thousand ulps of the resulting pmf. `n`, `1 - p` and
the products are therefore carried as double-doubles via `bd0_dd`. R's `dbinom`
has the same limitation.
Also guards `bd0`'s direct branch, where `(x / np).ln()` silently loses the whole
term when the ratio leaves the normal range - `bd0(1e-300, 5e99)` gave `-inf`
where the answer is `+5e99`, which propagated into `beta_reg` as NaN for extreme
parameter ratios, and as a silent `1.0` where the true value was ~1e-30118.
This was referenced Jul 26, 2026
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Densities and incomplete-function prefixes in the gamma/beta family are written as
exp(a·ln x − x − ln_gamma(a))and friends. Every term grows likea·ln awhile the sum stays O(ln a), so at a = 1e4 each term is ~1e5 and the ~1e-11 absolute error of the sum becomes the accuracy ceiling of everything downstream — this is what limitedbeta_regto ~6e-11 even where its recurrence converged, and putBinomial::pmfat ~7e-11 relative.Adds the two kernels from Loader, Fast and Accurate Computation of Binomial Probabilities (2000):
stirling_delta(z)— the Stirling-series remainder ofln_gamma, O(1/12z). Table-free: the recurrenceδ(z) = δ(z+1) + (z+½)·ln1p(1/z) − 1lifts any argument to where a 6-term series is good to 1.4e-18, so there is no interpolation table and no recursive dependency onln_gamma. ≤1.2e-15 absolute, which is what matters since every caller adds it to an exponent.bd0(x, np)=x·ln(x/np) − (x − np), with a double root atx = np. Evaluated by a series in(x−np)/(x+np)so the cancellation happens analytically: ~1e-16 relative uniformly through the root.Rewriting the prefixes (
gamma_lr/gamma_ur,beta_reg,Binomial::pmf,Poisson::pmf) in these terms keeps every piece O(1):beta_reg, a = b = 1e4beta_reg, a = b = 1e6gamma_ur, a ≥ 16Binomial::pmfPoisson::pmfTwo details reviewers should look at:
FisherSnedecor(1,1).cdfgot worse (the recurrence costs 16 steps at a = 0.5).bd0is sensitive enough that the rounding ofn·pitself matters (∂bd0/∂np = 1 − x/np: half an ulp ofnpat 6e5 moves the result 2.5e-13 ≈ 1000 ulp of the pmf), son,1−p, and the products are carried as double-doubles viabd0_dd. R'sdbinomhas this same limitation.Also guards
bd0's direct branch where(x/np).ln()under/overflows —bd0(1e-300, 5e99)gave−infwhere the answer is+5e99, which surfaced as NaN (or a silent 1.0 where the truth is ~1e-30118) inbeta_regat extreme parameter ratios.This is the same log-space reasoning as #381, applied to the rest of the family.
🤖 Generated with Claude Code