feat: add the inverse Gaussian (Wald) distribution#423
Open
agene0001 wants to merge 9 commits into
Open
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.
Adds `ln_cdf`/`ln_sf` to `ContinuousCDF` and `DiscreteCDF` with default implementations of `self.cdf(x).ln()`, plus tail-accurate overrides for the twelve distributions where the machinery exists. This is R's `log.p = TRUE`: every tail probability below ~1e-308 is currently irrecoverable - `Normal::sf(39)` is 0, so its log is -inf where the true value is -765.08 - and p-values at that scale are routine after multiple-testing correction (the use case in statrs-dev#338). The log is in most cases already computed internally and then exponentiated, so the overrides return the exponent instead of destroying it: * `erf::ln_erfc` (new): reuses the erfc interval machinery below 110, continues with the asymptotic series past it; finite until `x^2` overflows near 1.3e154. 0.34 median / 1.71 max ulp against mpmath over the full range. Backs Normal and LogNormal, including the argument-rounding compensation applied additively. * `gamma::ln_gamma_lr` / `ln_gamma_ur` (new, + checked variants): the prefix stays in log domain and the continued fraction / series are shared with the linear versions. Backs Gamma, ChiSquared, Erlang, Poisson. * `beta::ln_beta_reg` (new, + checked): likewise. Backs Beta and Binomial. * exact closed forms for Exp (`-rate x`), Weibull (`-(x/scale)^shape`), Pareto (`shape ln(scale/x)`), and Geometric (`x ln1p(-p)`). Validated against mpmath at 60 significant digits across all twelve distributions: the log-domain values carry the probability's *relative* accuracy at ~1e-15 uniformly, including regimes where cdf/sf are hard zeros (e.g. `Binomial(0.5, 1e4).ln_cdf(1000) = -3684.84...` where `cdf` is 0). The consistency test (`ln_cdf` vs `cdf().ln()` where both are representable) caught a pre-existing bug as a bonus: `Exp::cdf` used `1 - exp(-rate x)`, which cancels catastrophically for small `rate * x` (~7 digits lost at 1e-10). It now uses `-expm1`, matching what Weibull already did. StudentsT and FisherSnedecor keep the default implementation for now; their piecewise cdfs need a dedicated treatment. Closes statrs-dev#338
Closes statrs-dev#371. Implements the full trait surface: Continuous, ContinuousCDF (with the tail-accurate ln_cdf/ln_sf from the parent commit), Min, Max, Mode, Distribution's moments, and rand sampling via Michael-Schucany-Haas. The cdf is conventionally written Phi(sqrt(l/x) (x/m - 1)) + e^(2l/m) Phi(-sqrt(l/x) (x/m + 1)) which cannot be evaluated as written: e^(2l/m) overflows once l/m > ~355, even though the product it sits in is a probability bounded by 1. Both terms are therefore expressed through erfcx, using the identity 2l/m - u_plus^2 == -u_minus^2 (verified exactly in mpmath), so the large factor cancels symbolically and is never formed: cdf = e^(-u_minus^2) [erfcx(-u_minus) + erfcx(u_plus)] / 2 sf = e^(-u_minus^2) [erfcx( u_minus) - erfcx(u_plus)] / 2 whichever of the two is the smaller tail; the other is its complement. The e^(-u_minus^2) factor is Dekker-compensated for the rounding of u_minus^2, and the logs add the exponent exactly rather than taking ln of a product. This also adds erfcx and ln_erfcx to function::erf, which is what makes ratios of erfc values well conditioned (addresses statrs-dev#202). ln_erfcx is finite across the whole f64 range, where erfc underflows above ~27 and erfcx overflows below ~-27. Accuracy against mpmath at 60 digits, over 1800 points spanning six parameter sets and ten decades each. Wherever the probability is a representable nonzero f64: pdf 2.54 ulp median, 583 ulp max cdf 0.11 ulp median, 1179 ulp max sf 0.12 ulp median, 78686 ulp max (~9e-12 relative) ln_cdf 1.1e-23 median, 2.4e-13 max absolute, i.e. relative on ln_sf 9.4e-20 median, 1.0e-11 max the probability The linear-domain maxima are exp/cancellation limits rather than defects in these expressions: a relative error eps in an exponent becomes |exponent| * eps in the result, and sf differences two erfcx values whose ratio tends to (x/m - 1)/(x/m + 1). Past the underflow point the logs remain within 1-2 ulp of the returned value, where at ln = -3.4e8 one ulp is itself 6e-8. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
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.
Closes #371. Also adds
erfcx/ln_erfcx, which addresses #202.What's here
The full trait surface —
Continuous,ContinuousCDF(including tail-accurateln_cdf/ln_sf),Min,Max,Mode,Distribution's moments, andrandsampling via Michael-Schucany-Haas.
Why the cdf isn't written the textbook way
The conventional form is
which cannot be evaluated as written:
e^(2l/m)overflows oncel/m > ~355,even though the product it sits in is a probability bounded by 1. Both terms go
through
erfcxinstead, using(verified exactly in mpmath across six parameter sets) so the large factor
cancels symbolically and is never formed:
whichever of the two is the smaller tail; the other is returned as its
complement. The
e^(-u_minus^2)factor is Dekker-compensated for the roundingof
u_minus^2, and the log variants add that exponent exactly rather thantaking
lnof a product - an earlier draft differenced two logs of magnitude5e4to obtain2e-5, turning 0.7 ulp into4e-7of absolute error.erfcx(x) = e^(x^2) erfc(x)is also what makes ratios oferfcvalues wellconditioned generally (#202).
ln_erfcxis finite across the wholef64range, where
erfcunderflows above ~27 anderfcxoverflows below ~-27.Accuracy
Against mpmath at 60 digits over 1800 points - six parameter sets, ten decades
each. Restricted to points where the probability is a representable nonzero
f64, since that is the range a caller can observe a value in:pdfcdfsfln_cdfln_sfThe log rows are absolute error, which for a log is the relative error of the
probability it represents - so the implied probabilities are good to 12 and 11
significant digits respectively.
The linear-domain maxima are
expand cancellation limits, not defects inthese expressions: a relative error
epsin an exponent becomes|exponent| * epsin the result, andsfdifferences twoerfcxvalues whoseratio tends to
(x/m - 1)/(x/m + 1). Both are documented on the methods with apointer to the log variants.
Past the underflow point the logs stay within 1-2 ulp of the returned value.
Worth noting how that reads in absolute terms: the worst absolute error over the
unrestricted sweep is
1.2e-7, but it occurs atln cdf = -3.4e8, where oneulp is itself
6e-8. Neither metric is meaningful alone - ulp-of-the-logexplodes where the log approaches zero, and absolute error explodes where the
log is huge - so both are reported.
Verification
cargo fmt --check,cargo clippy --all-features --all-targetsclean,cargo hack check --feature-powerset --no-dev-depsclean,cargo hack check --rust-version --lockedclean againstCargo.lock.MSRV, and the full testsuite across
--all-features, default,--no-default-features, and--no-default-features --features rand.🤖 Generated with Claude Code