feat: add the inverse Gaussian (Wald) distribution - #423
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 is not a reproduction of the full 1800-point sweep from the PR. Representative values
Maximum observed errorRelative error is used for
Tail accuracy
The ordinary-domain results are effectively identical to SciPy. The significant improvement is in the deep right tail, especially PerformanceTimes are medians from repeated runs. SciPy vector timings are per element of a large array; scalar SciPy includes Python call overhead. mpmath uses 100-digit arithmetic and is intended only as a reference.
Compared with vectorized SciPy, the PR was approximately 1.4x–3.5x faster per element on this machine. Scalar SciPy is not a direct implementation-level comparison because it is dominated by Python API overhead. |
| /// Used instead of `f64::mul_add`, which falls back to a slow software FMA on | ||
| /// targets without the hardware instruction (e.g. baseline x86-64). | ||
| #[inline] | ||
| pub(crate) fn dekker_product_err(a: f64, b: f64, p: f64) -> f64 { |
There was a problem hiding this comment.
Why is the Dekker path unconditional? On AArch64/M1, a.mul_add(b, -p) compiles to a single fnmsub, while this implementation expands to roughly 30 FP instructions. The software fallback concern applies to targets without hardware FMA, such as generic x86-64, but this is a regression on FMA-capable targets.
Rust already selects hardware FMA based on target-cpu / target-feature, and the mul_add documentation explicitly describes this performance distinction. Could we use mul_add where hardware FMA is available and keep Dekker only as the fallback, with benchmarks for both paths?
| /// `-u^2` exactly, as the leading term of a log-domain result. | ||
| fn neg_square_exact(u: f64) -> f64 { | ||
| let sq = u * u; | ||
| -sq - crate::prec::dekker_product_err(u, u, sq) |
There was a problem hiding this comment.
This does not preserve the product residual. Since sq is already the rounded u * u, evaluating -sq - err rounds back to -sq before the caller adds the logarithmic term. For ordinary finite inputs this function is therefore effectively identical to -(u * u), so the log variants do not actually “add the exponent exactly”. The residual needs to be carried into the final sum, for example as (-sq + log_term) - err, or returned as a (hi, lo) pair.
Both from @day01's review on statrs-dev#423. `dekker_product_err` was applying Veltkamp's split unconditionally. It now dispatches: `mul_add` where the target has a hardware FMA, the split where it does not. Measured on an M1, 20M calls each: 0.83 ns against 2.48 ns, a 3.0x difference, with identical results. The reason for not using `mul_add` unconditionally still holds -- without hardware FMA it calls a software routine much slower than the split -- so the condition is the target, not the call site. Baseline x86-64 keeps the split; AArch64 has FMA in the base ISA. Both paths stay compiled everywhere and a test asserts they agree bit for bit over 200k random pairs plus the extremes, including the arguments past 1e300 where the split needs rescaling and the fused form does not. Otherwise the path not selected on the build host would never be exercised. The second point was right too, and my doc claim that the log variants "add the exponent exactly" was false. `neg_square_exact` computed `-sq - err` where `err` is at most half an ulp of `sq`, so the subtraction rounded straight back to `-sq` and the compensation did nothing. The suggested `(-sq + log) - err` does not fix it either -- the result still has magnitude `sq`, so `err` is still below half an ulp of it, and on every value I checked that ordering is bit-identical to the original. What does work is the other suggestion, carrying a pair: `neg_square_parts` returns `(-sq, -err)` unevaluated and `add_log_term` folds `lo` in through a two-sum with the log term. Over the 1800-point sweep this changes 37 of 360 ln_cdf points and halves the worst case, 2.4e-13 to 1.2e-13; ln_sf is unchanged at 1.0e-11. Docs updated to the measured figure. Also fixes a rustdoc `-D warnings` failure on this branch stack, where a public doc comment linked to the private `ln_beta_prefix`. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
| // The textbook root `mu + (mu / 2 lambda) (mnu - sqrt(mnu (4 lambda + | ||
| // mnu)))` cancels catastrophically for `mnu >> lambda`; rationalizing | ||
| // gives this equivalent, subtraction-free form. | ||
| let y = self.mu - (2.0 * self.mu * mnu) / (mnu + (mnu * (4.0 * self.lambda + mnu)).sqrt()); |
There was a problem hiding this comment.
z can be exactly zero because the RNG produces a finite set of f64 values. Substituting z = 0 gives nu = 0, mnu = 0, and:
y = mu - 0 / (0 + sqrt(0)) = mu - 0 / 0 = NaN.
| /// μ [sqrt(1 + 9μ^2 / (4λ^2)) - 3μ / (2λ)] | ||
| /// ``` | ||
| fn mode(&self) -> Option<f64> { | ||
| let r = self.mu / self.lambda; |
There was a problem hiding this comment.
| let r = self.mu / self.lambda; | |
| let ratio = self.lambda / self.mu; |
| /// ``` | ||
| fn mode(&self) -> Option<f64> { | ||
| let r = self.mu / self.lambda; | ||
| Some(self.mu * ((1.0 + 2.25 * r * r).sqrt() - 1.5 * r)) |
There was a problem hiding this comment.
| Some(self.mu * ((1.0 + 2.25 * r * r).sqrt() - 1.5 * r)) | |
| Some(self.lambda / ((ratio * ratio + 2.25).sqrt() + 1.5)) |
There was a problem hiding this comment.
This is algebraically equivalent to the documented mode formula, but avoids catastrophic cancellation.
example mu = 1e8 and lambda = 1:
current:
r = mu / lambda = 1e8
mode = 1e8 * (sqrt(1 + 2.25 * 1e16) - 1.5e8)
= 1e8 * (1.5e8 - 1.5e8)
= 0
suggests:
r = lambda / mu = 1e-8
mode = 1 / (sqrt(1e-16 + 2.25) + 1.5)
≈ 1 / 3
| /// a `-C target-cpu` that implies it -- while AArch64 has it in the base ISA. | ||
| #[inline] | ||
| pub(crate) fn dekker_product_err(a: f64, b: f64, p: f64) -> f64 { | ||
| #[cfg(any(target_arch = "aarch64", target_feature = "fma"))] |
There was a problem hiding this comment.
| #[cfg(any(target_arch = "aarch64", target_feature = "fma"))] | |
| #[cfg(any( | |
| target_feature = "fma", | |
| all(target_arch = "aarch64", target_feature = "neon") | |
| ))] |
aarch as softfloat has software fma n there dekker is good choice
There was a problem hiding this comment.
Shouldn't the FMA path be used on every target where mul_add lowers to a hardware instruction? It should generally be faster than the Dekker split there, not only on x86 with fma and AArch64.
day01
left a comment
day01
left a comment
There was a problem hiding this comment.
if we re adding fma sup, it shouldnt be only for aarch, it is wider
|
i have a question @YeungOnion , if @agene0001 is changing erfcx, so maybe better to user erfcx as implementation of Cody? I made it times ago: https://github.com/day01/quant-opts/blob/ed23cd715da46da3b8cd52de09904e95b5b1a967/src/lets_be_rational/cody/optimized.rs#L119 imo it is faster and accuracy is better, but maybe i am biased. |
All three from @day01's review on statrs-dev#423. **Mode.** `sqrt(1 + 9u^2/4L^2) - 3u/2L` differences two terms that converge as `mu/lambda` grows, so it cancels: at `mu = 1e8, lambda = 1` both are `1.5e8` to working precision and `mode()` returned `0` instead of `1/3`. Now evaluated as `lambda / (sqrt((lambda/mu)^2 + 9/4) + 3/2)`, which is the same expression multiplied through by its conjugate. That also keeps the squared term from overflowing, so `mu = 1e300` works where it previously could not. Worth noting how the references for the new test were produced: evaluating the textbook form in mpmath at 50 digits returns zero for `mu = 1e150` too. It needs about 400 to give an answer at all. The reference values come from the conjugate form. **The sampler could return NaN.** `z` really can be exactly zero -- the ziggurat forms `u = 2f - 1` from a 53-bit `f`, and `f == 0.5` is one of the values it takes, so this has probability `2^-53` per draw rather than being impossible. The rationalized root then divided by `mnu + sqrt(mnu (4 lambda + mnu))`, which is `0/0` there, and the NaN propagated through both branches of the conjugate selection. Fixed by pulling `sqrt(mnu)` out of the radical, which cancels it against the numerator and leaves a denominator bounded below by `2 sqrt(lambda)`. The result is finite everywhere and tends to `mu`, as it should for `nu = 0`. It also stops `mnu (4 lambda + mnu)` overflowing for large `mnu`. Split out as `msh_smaller_root` so the boundary is testable without an RNG, since no practical number of draws would reach it. **FMA detection was too narrow**, on both counts raised. AArch64 was treated as unconditional, but a softfloat AArch64 target has no NEON and so no hardware FMA, where the split is the right choice; and only x86 was considered besides. Detection moves to build.rs, which sets `cfg(statrs_hardware_fma)` from the target's actual features, covering AArch64, x86, ARMv7, RISC-V, PowerPC, LoongArch, MIPS and WASM. Each arm tests for the instruction rather than assuming the architecture implies it. Reverting either numerical fix fails its test. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
|
All three addressed. Both numerical findings were real — thanks, these were good catches. ModeConfirmed exactly as you described: One thing worth adding, which I hit while producing test references: evaluating the textbook form in mpmath at 50 digits also returns zero for The NaNAlso real, and reachable rather than theoretical. The ziggurat builds Rather than special-casing it, pulling The denominator is now bounded below by FMA detectionYou were right on both points, and the second one is the one I got wrong more broadly: I had treated AArch64 as unconditional and considered nothing outside x86. Detection now lives in
Each arm tests for the instruction rather than assuming the architecture implies it, so softfloat AArch64 correctly takes the split. The equivalence test between the two paths is unchanged and still runs both everywhere. On the Cody
|
Found by benchmarking against @day01's Cody erfcx, which was more accurate here by two orders of magnitude and prompted the question of why. The branch computed `exp(x * x) * erfc(x)`. `x * x` carries up to half an ulp and exponentiating multiplies that error by `x^2`, so the result drifted by roughly `x^2 / 2` ulp -- about 340 by `x = -26`, where the branch bottoms out before `exp` overflows. The comment claimed `exp(x^2)` was "harmless up to x^2 = 709", which is true of overflow and not of accuracy. Recovering the residual with the Dekker split already used elsewhere in this file, and folding it back as `exp(e) ~ 1 + e`, removes the amplification. Over 7604 points against mpmath at 60 digits, for `x < -1`: median 31.04 -> 0.38 ulp max 463.08 -> 1.77 ulp The other branches are untouched. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
|
@day01 I ran the comparison. You were right on both counts, and it found a real bug in mine — thanks for pushing on this. Setup: 7604 points against mpmath at 60 digits, covering Accuracy, as it stood
Yours was ~140x better on the worst case for WhyMy negative branch was Your
Overall, statrs is now 0.36 median / 2.08 max against Cody's 0.40 / 4.21, and better in every bucket. SpeedYours is still clearly faster:
The compensation cost me some of that (7.64 → 10.26 ns on the mixed path); the positive path is unchanged. SoThe honest summary is that it's now an accuracy-versus-speed trade rather than a clear win either way: statrs is ~2x more accurate on worst case, Cody is ~2.3x faster. Which matters more is a call for @YeungOnion — for a stats library I'd lean accuracy, but I'm obviously not neutral, and 2x throughput on a hot path is not nothing. Either way the bug your question exposed is fixed independently of that decision. One methodological note in case it's useful: my first reference values for the regression test were wrong, because I generated them from decimal literals. |
2 participants
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