fix: compensate the argument rounding in Normal's cdf/sf tails#419
Draft
agene0001 wants to merge 7 commits into
Draft
fix: compensate the argument rounding in Normal's cdf/sf tails#419agene0001 wants to merge 7 commits into
agene0001 wants to merge 7 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.
`erf_impl` reconstructs `erfc(z)` for `z >= 0.5` as
`exp(-z^2)/z * (b + r(z))`, where `r` is a minimax rational fit and `b` is a
per-interval constant. All 13 of those constants carried only 10 significant
decimal digits:
0.3440242112 // interval [0.5, 0.75)
0.419990927 // interval [0.75, 1.25)
Solving for the `b` each interval's fit implies, at 60 digits, gives a value
that is constant to ~19 digits across the whole interval - so the rational fits
are good to ~1e-19 and the truncated constants were the only error source. The
recovered values turn out to be exactly `f32`-representable, i.e. Boost declared
them as `float` literals and the port that statrs inherits transcribed the
printed decimal:
0.3440242112 -> 0.3440242111682891845703125
0.419990927 -> 0.4199909269809722900390625
Also compensates the squaring: `z * z` is rounded before `exp` sees it, which
costs roughly `z^2` ulps (~460 by z = 27). A Dekker product recovers the
rounding error of the square exactly and folds it back in as
`exp(-z^2) = exp(-sq) * (1 - err)`. (Dekker rather than `f64::mul_add`, which
falls back to a slow software FMA on targets without the instruction.)
Measured against mpmath at 45 dps over a 700-point sweep plus every interval
boundary:
before after
erf ~1e-10 rel 0.24 median / 1.1 max ulp
erfc ~1e-10 rel 0.52 median / 2.8 max ulp
Normal::cdf ~5e-11 rel 0.32 median
This propagated to everything built on `erfc` - `Normal` cdf/sf, `LogNormal`,
`Levy`. The test suite was masking it with `epsilon = 1e-9`/`1e-11`
tolerances, now tightened to 4 ulp.
Six reference literals in the tests were themselves wrong, having been fitted to
the old output; they are replaced with mpmath values. Three are `LogNormal::sf`
expectations (off by up to 5.1e-12) and three are `erfc_inv` expectations - one
of which, `erfc_inv(1e-10)`, was off by 8.8e-9 and had needed
`epsilon = 1e-7` to pass. `erf_inv`/`erfc_inv` themselves were already good to
<2 ulp; only the expectations were wrong.
`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.
`Normal::cdf` is `0.5 * erfc((mean - x) / (sigma * sqrt 2))`, and `erfc`
amplifies a relative error in its argument by `2 t^2`. The subtraction, the
division, and the representation of `sigma * sqrt(2)` each contribute a rounding,
so the tails carried ~40-80 ulp even once `erfc` itself was correct.
Recovers all three residuals exactly - a two-sum for the difference, a Dekker
product for the division, and a double-double `sqrt(2)` - and applies the
first-order correction `erfc(t + d) ~= erfc(t) + erfc'(t) d`. `erfc'` needs
`exp(-t^2)`, which is the only added cost, so it is applied only for
`1.5 < |t| < 30`: below that the amplification is negligible and `erfc` is
absolutely well conditioned, above it `erfc` has already saturated to exactly 0
or 2.
Measured against mpmath at 45 dps over 601 points per parameter set:
N(0, 1) 80.8 max ulp -> 4.0
N(0.3, 1.7) ~80 -> 7.5
N(-2.5, 0.013) -> 5.6
N(1e4, 250) -> 3.3
Centre calls are unchanged at 8.7 ns; tail calls cost 15.8 ns for the extra
`exp`.
The same range gate also keeps the Dekker splits away from arguments they cannot
handle - they multiply by `2^27 + 1`, which overflows past ~1.3e300 and would
turn the result into NaN - and infinite `x` or an overflowing `sigma * sqrt(2)`
are handled explicitly, since `inf / inf` is NaN where both limits are well
defined. `N(0, f64::MIN_POSITIVE).cdf(-1)` and `N(0, f64::MAX).cdf(-inf)` were
NaN; a test now asserts cdf/sf stay in `[0, 1]` across seven parameter sets and
nine inputs including `+-f64::MAX`.
Depends on the erf constant fix (merged in), since compensating the argument of a
1e-10-accurate `erfc` would be pointless.
agene0001
marked this pull request as draft
July 27, 2026 08:32
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.
Normal::cdfis0.5·erfc((mean − x)/(σ·√2)), anderfcamplifies a relative error in its argument by2t². The subtraction, the division, and the representation ofσ·√2each contribute a rounding, so the tails carried ~40–80 ulp even onceerfcitself was accurate.half_erfcnow recovers all three residuals exactly — a two-sum for the difference, a Dekker product for the division, and a double-double√2— and applies the first-order correctionerfc(t + δ) ≈ erfc(t) + erfc′(t)·δ. The correction needs one extraexp, so it is applied only for1.5 < |t| < 30: below that the amplification is negligible anderfcis absolutely well-conditioned; above iterfchas saturated to exactly 0 or 2.Measured against mpmath at 45 dps, 601 points per parameter set:
Centre calls unchanged (8.7 ns); tail calls 15.8 ns.
The same gate keeps the Dekker splits away from arguments they can't handle (they multiply by 2²⁷+1, which overflows past ~1.3e300 into NaN), and infinite
x/ overflowingσ·√2are handled explicitly sinceinf/infis NaN where both limits are well-defined.N(0, f64::MIN_POSITIVE).cdf(-1)andN(0, f64::MAX).cdf(-inf)were NaN; a test now asserts cdf/sf stay in[0, 1]across seven parameter sets × nine inputs including ±f64::MAX.🤖 Generated with Claude Code