diff --git a/src/factor/trunc.rs b/src/factor/trunc.rs index f2d04a3..3a4184d 100644 --- a/src/factor/trunc.rs +++ b/src/factor/trunc.rs @@ -2,6 +2,7 @@ use crate::{ N_INF, approx, cdf, factor::{Factor, VarId, VarStore}, gaussian::Gaussian, + sf, }; /// EP truncation factor on a diff variable. @@ -74,16 +75,29 @@ impl Factor for TruncFactor { /// P(diff > margin) for non-tie, P(|diff| < margin) for tie. /// -/// Clamped to a positive floor: for a near-certain outcome the tail rounds to -/// exactly 0.0, and the `erfc` approximation used by `cdf` carries ~1e-7 error -/// so it can even return slightly more than 1.0, making the difference -/// negative. Either would send `log_evidence` to `-inf` or NaN and poison the -/// sum across the whole history. +/// Both branches pick whichever tail keeps their terms *small*, because the +/// alternative is subtracting two numbers that both approach 1. That +/// subtraction is not a rounding detail: it loses every digit of an unlikely +/// outcome's evidence, and an unlikely outcome is precisely the one worth +/// scoring. `1 - cdf` returned exactly zero past ~8.3 sigma, where the true +/// probability is 1e-19; clamped, that reached `log_evidence` as -708 instead +/// of -43. +/// +/// The clamp remains as a guard rather than a workaround: `erfc` carries ~1e-7 +/// relative error, so a probability of exactly 1 can still come back a hair +/// above it, and `ln` of a negative would poison the sum for the whole history. fn cavity_evidence(diff: Gaussian, margin: f64, tie: bool) -> f64 { + let (mu, sigma) = (diff.mu(), diff.sigma()); + let raw = if tie { - cdf(margin, diff.mu(), diff.sigma()) - cdf(-margin, diff.mu(), diff.sigma()) + if mu < -margin { + // Both CDFs sit against 1 here; both survival terms are small. + sf(-margin, mu, sigma) - sf(margin, mu, sigma) + } else { + cdf(margin, mu, sigma) - cdf(-margin, mu, sigma) + } } else { - 1.0 - cdf(margin, diff.mu(), diff.sigma()) + sf(margin, mu, sigma) }; raw.clamp(f64::MIN_POSITIVE, 1.0) @@ -132,6 +146,52 @@ mod tests { assert_eq!(f.evidence_cached.unwrap(), first); } + /// The defect this guards: `1 - cdf` collapsed to zero for a surprising + /// result, the clamp turned that into `f64::MIN_POSITIVE`, and + /// `log_evidence` reported ln of *that* — about -708 whatever the truth + /// was. An upset is the observation a model-comparison score exists to + /// notice, so it was wrong exactly where it mattered. + #[test] + fn evidence_of_an_upset_is_not_flattened_to_the_clamp_floor() { + // diff ~ N(-9, 1) with margin 0: the favoured side lost by nine sigma. + let evidence = cavity_evidence(Gaussian::from_ms(-9.0, 1.0), 0.0, false); + + assert!( + evidence > f64::MIN_POSITIVE, + "evidence collapsed onto the clamp floor: {evidence}" + ); + // P(X > 0) for X ~ N(-9, 1) is the standard normal tail at 9 sigma. + assert!( + (evidence - 1.128_588e-19).abs() / 1.128_588e-19 < 1e-6, + "expected ~1.13e-19, got {evidence}" + ); + assert!( + (evidence.ln() + 43.628).abs() < 1e-2, + "log evidence {} should be about -43.6, not -708", + evidence.ln() + ); + } + + /// Evidence must stay finite and positive however extreme the mismatch, + /// since `log_evidence` sums across the whole history and one `-inf` or + /// `NaN` poisons all of it. + #[test] + fn evidence_stays_positive_and_finite_at_any_separation() { + for mu in [-300.0f64, -50.0, -9.0, 0.0, 9.0, 50.0, 300.0] { + for tie in [false, true] { + let e = cavity_evidence(Gaussian::from_ms(mu, 1.0), 1.0, tie); + assert!( + e.is_finite() && e > 0.0 && e <= 1.0, + "mu={mu} tie={tie}: evidence {e} is not a probability" + ); + assert!( + e.ln().is_finite(), + "mu={mu} tie={tie}: ln evidence is not finite" + ); + } + } + } + #[test] fn tie_evidence_uses_two_sided() { let mut vars = VarStore::new(); diff --git a/src/lib.rs b/src/lib.rs index 01a6a5a..40a8762 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -168,6 +168,21 @@ pub const ITERATIONS: usize = 30; pub const MAX_PREDICTED_TEAMS: usize = predict::MAX_TEAMS_FOR_DISTRIBUTION; const SQRT_TAU: f64 = 2.5066282746310002; +/// `1 / sqrt(pi)`, the leading factor of the `erfcx` continued fraction. +const FRAC_1_SQRT_PI: f64 = 0.564_189_583_547_756_3; +/// `sqrt(2 / pi)`, the numerator of the inverse Mills ratio in scaled form. +const SQRT_2_OVER_PI: f64 = 0.797_884_560_802_865_4; +/// How many window widths into the tail before a tie window is treated as a +/// half-line. Beyond this the truncated mass is concentrated within `1/alpha` +/// of the near edge, so the far edge contributes nothing measurable. +const HALF_LINE_WINDOW: f64 = 10.0; +/// Where `v - alpha` switches from subtraction to its asymptotic series. +/// +/// The subtraction loses roughly `eps * alpha^2` of relative precision, and the +/// four-term series is good to ~1e-10 by here, so the two are at their closest +/// agreement around this point. Below it the subtraction is exact; above it the +/// series is. +const ASYMPTOTIC_MILLS_ALPHA: f64 = 100.0; pub const N01: Gaussian = Gaussian::from_ms(0.0, 1.0); pub const N00: Gaussian = Gaussian::from_ms(0.0, 0.0); @@ -259,6 +274,50 @@ pub(crate) fn cdf(x: f64, mu: f64, sigma: f64) -> f64 { 0.5 * erfc(z) } +/// `P(X > x)` for `X ~ N(mu, sigma^2)`. +/// +/// The survival function, computed directly rather than as `1 - cdf(..)`. +/// +/// The two are algebraically identical and numerically are not. `cdf` returns +/// a value approaching 1 for an upper tail, so subtracting it from 1 cancels +/// away every significant digit the tail had: measured against this function, +/// `1 - cdf` carries 7% error by four sigma past the mean and returns exactly +/// zero beyond about 8.3 sigma — where the true value is still 1e-19 and +/// perfectly representable. `erfc` itself holds ~1e-7 *relative* accuracy down +/// to 1e-296, so the precision is there to keep; only the subtraction threw it +/// away. +/// +/// This matters most where evidence is smallest, which is exactly where an +/// upset makes it interesting: `ln` of a clamped zero is -708 regardless of +/// whether the truth was -43 or -600. +pub(crate) fn sf(x: f64, mu: f64, sigma: f64) -> f64 { + 0.5 * erfc((x - mu) / (sigma * SQRT_2)) +} + +/// `e^(x^2) * erfc(x)`, the scaled complementary error function, for `x >= 0`. +/// +/// Exists so the exponential factor common to a Gaussian density and its tail +/// integral can be cancelled *analytically* instead of being computed twice +/// and divided. Both underflow to zero past about 26 sigma, and their ratio is +/// then `0/0` — finite in the limit, `NaN` in floating point. +fn erfcx(x: f64) -> f64 { + if x < 2.0 { + // Below the crossover neither factor is extreme: erfc is O(1) and + // exp(x^2) is at most e^4, so the direct product is exact enough and + // cheaper than the continued fraction. + (x * x).exp() * erfc(x) + } else { + // erfcx(x) = 1/sqrt(pi) * 1/(x + (1/2)/(x + 1/(x + (3/2)/(x + ...)))), + // evaluated by backward recurrence. Converges quickly for x >= 2 and, + // unlike the product form, never touches an exponential. + let mut f = 0.0; + for n in (1..=60u32).rev() { + f = (f64::from(n) * 0.5) / (x + f); + } + FRAC_1_SQRT_PI / (x + f) + } +} + fn pdf(x: f64, mu: f64, sigma: f64) -> f64 { let normalizer = (SQRT_TAU * sigma).powi(-1); let functional = (-((x - mu).powi(2)) / (2.0 * sigma.powi(2))).exp(); @@ -266,25 +325,100 @@ fn pdf(x: f64, mu: f64, sigma: f64) -> f64 { normalizer * functional } +/// Truncated-Gaussian correction terms `(v, w)`. +/// +/// `v` shifts the mean and `w` shrinks the variance. Both are ratios whose +/// numerator and denominator underflow together in the tails, so both are +/// computed in scaled form there: the shared `exp(-alpha^2 / 2)` is cancelled +/// analytically rather than evaluated and divided out. Without that, a +/// truncation point beyond about 39 sigma produced `0 / 0` and put `NaN` +/// straight into the posterior. +/// Truncation terms for a boundary `alpha` standard deviations into the upper +/// tail, from the asymptotic expansion of the inverse Mills ratio. +/// +/// `v` tends to `alpha` out here, so the gap between them cannot be obtained by +/// subtracting one from the other — the series computes the gap directly, and +/// `w = v * gap` then never forms the difference of two large near-equal +/// numbers. A far-tail *window* behaves like a half-line once it is more than a +/// few multiples of its own width from the mean, so the tie branch shares this. +fn half_line_truncation(alpha: f64) -> (f64, f64) { + let inv = alpha.recip(); + let inv_sq = inv * inv; + let gap = inv * (1.0 - inv_sq * (2.0 - inv_sq * (10.0 - 74.0 * inv_sq))); + let v = alpha + gap; + + (v, v * gap) +} + fn v_w(mu: f64, sigma: f64, margin: f64, tie: bool) -> (f64, f64) { if !tie { let alpha = (margin - mu) / sigma; - let v = pdf(-alpha, 0.0, 1.0) / cdf(-alpha, 0.0, 1.0); - let w = v * (v + (-alpha)); + // v is the inverse Mills ratio, phi(alpha) / Phi(-alpha), and w needs + // the gap `v - alpha` as well as v itself. Far into the tail v tends to + // alpha, so that gap is a subtraction of two nearly equal numbers and + // loses every digit it has: at alpha = 1e6 it drove w above 1 and made + // `sqrt(1 - w)` NaN. Past the crossover the gap comes from its + // asymptotic series instead, which has no subtraction in it. + if alpha >= ASYMPTOTIC_MILLS_ALPHA { + return half_line_truncation(alpha); + } - (v, w) + let (v, gap) = if alpha > 0.0 { + // Both terms carry exp(-alpha^2 / 2); in scaled form it cancels + // and the result stays exact however far into the tail alpha sits. + let v = SQRT_2_OVER_PI / erfcx(alpha / SQRT_2); + (v, v - alpha) + } else { + // Phi(-alpha) >= 1/2 here, so the direct ratio loses nothing. + let v = pdf(-alpha, 0.0, 1.0) / cdf(-alpha, 0.0, 1.0); + (v, v - alpha) + }; + + (v, v * gap) } else { + // v is odd in mu and w is even, so fold to mu <= 0. Both truncation + // points then sit in the upper tail, where the scaled form applies. + let flipped = mu > 0.0; + let mu = if flipped { -mu } else { mu }; + let alpha = (-margin - mu) / sigma; let beta = (margin - mu) / sigma; - let v = (pdf(alpha, 0.0, 1.0) - pdf(beta, 0.0, 1.0)) - / (cdf(beta, 0.0, 1.0) - cdf(alpha, 0.0, 1.0)); - let u = (alpha * pdf(alpha, 0.0, 1.0) - beta * pdf(beta, 0.0, 1.0)) - / (cdf(beta, 0.0, 1.0) - cdf(alpha, 0.0, 1.0)); + // `w` comes out of `v * v - u`, and both terms grow as alpha^2 while + // their difference stays O(1) — at alpha = 1e9 that subtraction had no + // digits left and returned w = -128, making `sqrt(1 - w)` nonsense. + // Once the window sits many of its own widths into the tail it is + // indistinguishable from a half-line, so the asymptotic covers it with + // no subtraction at all. + if alpha >= ASYMPTOTIC_MILLS_ALPHA && alpha * (beta - alpha) >= HALF_LINE_WINDOW { + let (v, w) = half_line_truncation(alpha); + return (if flipped { -v } else { v }, w); + } + + let (v, u) = if alpha > 0.0 { + // beta > alpha > 0, so this ratio of exponentials is at most 1 and + // cannot overflow. + let scale = (0.5 * (alpha * alpha - beta * beta)).exp(); + let denominator = 0.5 * (erfcx(alpha / SQRT_2) - scale * erfcx(beta / SQRT_2)); + + ( + (1.0 - scale) / SQRT_TAU / denominator, + (alpha - beta * scale) / SQRT_TAU / denominator, + ) + } else { + // The interval straddles the mean, so nothing here is small. + let denominator = cdf(beta, 0.0, 1.0) - cdf(alpha, 0.0, 1.0); + + ( + (pdf(alpha, 0.0, 1.0) - pdf(beta, 0.0, 1.0)) / denominator, + (alpha * pdf(alpha, 0.0, 1.0) - beta * pdf(beta, 0.0, 1.0)) / denominator, + ) + }; + let w = -(u - v.powi(2)); - (v, w) + (if flipped { -v } else { v }, w) } } @@ -467,6 +601,151 @@ mod tests { assert_eq!(sort_time(&[0i64, 1, 2, 0], true), vec![2, 1, 0, 3]); } + /// Upper-tail values of the standard normal, from published tables. The + /// point is not the digits — `erfc` only carries ~1e-7 relative — but that + /// a number comes back at all: `1 - cdf` returned exactly zero for every + /// one of these. + #[test] + fn survival_function_survives_the_far_tail() { + for (z, expected) in [ + (9.0f64, 1.128_588e-19), + (12.0, 1.776_482e-33), + (20.0, 2.753_624e-89), + (37.0, 5.725_571e-300), + ] { + let got = sf(z, 0.0, 1.0); + assert!(got > 0.0, "sf({z}) collapsed to zero"); + assert!( + (got - expected).abs() / expected < 1e-6, + "sf({z}) = {got}, expected ~{expected}" + ); + assert_eq!( + 1.0 - cdf(z, 0.0, 1.0), + 0.0, + "the naive form should still be zero here" + ); + } + } + + /// Where no cancellation happens the two forms must agree exactly enough + /// that nothing else in the crate shifts. + #[test] + fn survival_function_matches_the_naive_form_where_that_form_works() { + for z in [-4.0f64, -1.0, 0.0, 0.5, 1.0, 2.0, 3.0, 4.0] { + let naive = 1.0 - cdf(z, 0.0, 1.0); + let direct = sf(z, 0.0, 1.0); + // Bounded by `erfc`'s own ~1e-7 relative error, not by the + // subtraction: the two forms evaluate `erfc` at different points + // and the approximation is not exactly antisymmetric. + assert!( + (naive - direct).abs() < 1e-6, + "z={z}: naive {naive} vs direct {direct}" + ); + } + } + + #[test] + fn survival_and_cdf_partition_the_mass() { + for z in [-3.0f64, -0.5, 0.0, 1.0, 2.5] { + let total = sf(z, 1.0, 2.0) + cdf(z, 1.0, 2.0); + // `erfc(z) + erfc(-z) == 2` only to the accuracy of the + // approximation, which is ~1e-7 relative. + assert!((total - 1.0).abs() < 1e-6, "z={z}: {total}"); + } + } + + /// `erfcx` switches formulation at x = 2; the two sides must meet. + #[test] + fn erfcx_is_continuous_across_its_crossover() { + for x in [1.90f64, 1.99, 1.999, 2.0, 2.001, 2.01, 2.10] { + let direct = (x * x).exp() * erfc(x); + let scaled = erfcx(x); + assert!( + (direct - scaled).abs() / scaled < 1e-6, + "x={x}: direct {direct} vs erfcx {scaled}" + ); + } + } + + /// The whole reason `erfcx` exists: it stays finite and O(1/x) exactly + /// where `exp(x^2)` overflows and `erfc(x)` underflows. + #[test] + fn erfcx_stays_finite_where_its_factors_do_not() { + for x in [27.0f64, 50.0, 1.0e3, 1.0e8] { + let scaled = erfcx(x); + assert!(scaled.is_finite() && scaled > 0.0, "erfcx({x}) = {scaled}"); + // Asymptotically erfcx(x) -> 1 / (x * sqrt(pi)). + let asymptote = 1.0 / (x * std::f64::consts::PI.sqrt()); + assert!( + (scaled - asymptote).abs() / asymptote < 1e-2, + "erfcx({x}) = {scaled} strays from its asymptote {asymptote}" + ); + assert!( + (x * x).exp().is_infinite(), + "x={x} should overflow the direct form" + ); + } + } + + /// Truncation must never produce a non-finite posterior. Before the scaled + /// formulation these returned NaN from `0 / 0` past about 39 sigma. + #[test] + fn truncation_stays_finite_arbitrarily_far_into_the_tail() { + for alpha in [0.0f64, 8.0, 38.0, 40.0, 100.0, 1.0e3, 1.0e6, 1.0e9, 1.0e15] { + for tie in [false, true] { + let (v, w) = v_w(-alpha, 1.0, if tie { 1.0 } else { 0.0 }, tie); + assert!(v.is_finite(), "alpha={alpha} tie={tie}: v = {v}"); + assert!(w.is_finite(), "alpha={alpha} tie={tie}: w = {w}"); + // sigma_trunc = sigma * sqrt(1 - w) must stay real. + assert!( + (0.0..=1.0).contains(&w), + "alpha={alpha} tie={tie}: w = {w} leaves sqrt(1 - w) imaginary" + ); + + let (mu_t, sigma_t) = trunc(-alpha, 1.0, if tie { 1.0 } else { 0.0 }, tie); + assert!( + mu_t.is_finite() && sigma_t.is_finite(), + "alpha={alpha} tie={tie}: trunc = ({mu_t}, {sigma_t})" + ); + } + } + } + + /// The Mills gap switches from subtraction to series at alpha = 100. Both + /// are supposed to be right there; if they disagree, the crossover is in + /// the wrong place. + #[test] + fn the_mills_gap_series_meets_the_scaled_form() { + for alpha in [50.0f64, 99.0, 100.0, 101.0, 200.0] { + let scaled = SQRT_2_OVER_PI / erfcx(alpha / SQRT_2) - alpha; + let inv = alpha.recip(); + let inv_sq = inv * inv; + let series = inv * (1.0 - inv_sq * (2.0 - inv_sq * (10.0 - 74.0 * inv_sq))); + assert!( + (scaled - series).abs() / series < 1e-9, + "alpha={alpha}: scaled {scaled} vs series {series}" + ); + } + } + + /// Folding the tie branch to `mu <= 0` is only valid if v is odd in mu and + /// w is even. Assert the symmetry the implementation relies on. + #[test] + fn tie_truncation_is_odd_in_v_and_even_in_w() { + for mu in [0.5f64, 3.0, 20.0, 40.0, 100.0, 1.0e3] { + let (v_pos, w_pos) = v_w(mu, 1.0, 1.0, true); + let (v_neg, w_neg) = v_w(-mu, 1.0, 1.0, true); + assert!( + (v_pos + v_neg).abs() < 1e-9, + "mu={mu}: v should be odd, got {v_pos} and {v_neg}" + ); + assert!( + (w_pos - w_neg).abs() < 1e-9, + "mu={mu}: w should be even, got {w_pos} and {w_neg}" + ); + } + } + #[test] fn test_quality() { let a = Gaussian::from_ms(25.0, 3.0);