fix: scale beta_reg's continued-fraction bound with its parameters#417
Open
agene0001 wants to merge 3 commits into
Open
fix: scale beta_reg's continued-fraction bound with its parameters#417agene0001 wants to merge 3 commits into
agene0001 wants to merge 3 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.
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.
The Lentz recurrence in
checked_beta_regwas capped at a fixed 140 iterations and silently returned whatever it had reached. The recurrence is slowest at the centre of the distribution, where the worst case overxgrows like5·min(a,b)^(1/3)(measured over a,b ∈ [1e2, 1e18] × 4001 values of x), so pastmin(a,b) ≈ 1.5e4the answer was simply wrong. Against the exact identityI_½(a,a) = ½:Binomial::new(0.5, 2e6).cdf(1e6)returned 0.4916 against a true 0.50028 — this capsStudentsT,FisherSnedecor, andBetacdfs too.The bound now scales as
8·min(a,b)^(1/3)with headroom, clamped to bound worst-case work at ~10 ms. The loop still exits on convergence, so ordinary calls are unchanged (beta_reg(2.5, 2.5, 0.5): 150 ns vs 145 ns).Plus four robustness fixes found by probing the parameter space rather than sweeping accuracy:
0.0 × h, which is NaN whenever the recurrence overflowed (beta_reg(1e300, 1e-300, 0.5)was NaN before, on both sides of the bound change);x == 0/x == 1return 0/1 directly instead of via the symmetry test, which mappedx == 0to 1.0 oncea + boverflowed;(a+1)/(a+b+2)is computed scaled whena + boverflows, otherwise it collapses to 0 and sends everyxdown the transformed branch;[0, 1](it's a probability;a = b = 1e20previously returned −2.56), falling back to the concentrated-limit step function if the truncated recurrence went non-finite.Regression tests use exact identities needing no reference data —
I_½(a,a) = ½andI_x(a,b) + I_{1−x}(b,a) = 1(both fail on the old fixed bound) — plus a 12×12×6 extreme-parameter grid asserting the result stays a finite probability.🤖 Generated with Claude Code