diff --git a/src/factor/margin.rs b/src/factor/margin.rs index 3721552..2d96900 100644 --- a/src/factor/margin.rs +++ b/src/factor/margin.rs @@ -2,7 +2,7 @@ use crate::{ N_INF, factor::{Factor, VarId, VarStore}, gaussian::Gaussian, - pdf, + ln_pdf, }; /// Gaussian observation factor on a diff variable. @@ -16,7 +16,7 @@ pub struct MarginFactor { pub m_obs: f64, pub sigma: f64, pub(crate) msg: Gaussian, - pub(crate) evidence_cached: Option, + pub(crate) log_evidence_cached: Option, } impl MarginFactor { @@ -28,7 +28,7 @@ impl MarginFactor { m_obs, sigma, msg: N_INF, - evidence_cached: None, + log_evidence_cached: None, } } } @@ -41,8 +41,8 @@ impl MarginFactor { let marginal = vars.get(self.diff); let cavity = marginal / self.msg; - if self.evidence_cached.is_none() { - self.evidence_cached = Some(cavity_evidence(cavity, self.m_obs, self.sigma)); + if self.log_evidence_cached.is_none() { + self.log_evidence_cached = Some(cavity_log_evidence(cavity, self.m_obs, self.sigma)); } let new_msg = Gaussian::from_ms(self.m_obs, self.sigma); @@ -61,17 +61,29 @@ impl Factor for MarginFactor { } fn log_evidence(&self, _vars: &VarStore) -> f64 { - self.evidence_cached.unwrap_or(1.0).ln() + self.log_evidence_cached.unwrap_or(0.0) } } -/// Density of the observed margin under the cavity, clamped to a positive -/// floor so a far-out observation cannot underflow to `0.0` and make -/// `log_evidence` `-inf`. -fn cavity_evidence(cavity: Gaussian, m_obs: f64, sigma: f64) -> f64 { +/// `ln` of the observed margin's density under the cavity. +/// +/// Computed in log space rather than as `pdf(..).ln()`. The density underflows +/// to zero past about 38 sigma of separation, and clamping that to +/// `f64::MIN_POSITIVE` reported -708 nats however far out the observation +/// actually was — 4292 nats adrift at 100 sigma, and unbounded beyond. A score +/// far from what the model expected is exactly the observation a log-evidence +/// figure exists to notice. +fn cavity_log_evidence(cavity: Gaussian, m_obs: f64, sigma: f64) -> f64 { let combined_sigma = (cavity.sigma().powi(2) + sigma.powi(2)).sqrt(); + let value = ln_pdf(m_obs, cavity.mu(), combined_sigma); - pdf(m_obs, cavity.mu(), combined_sigma).max(f64::MIN_POSITIVE) + // A degenerate cavity (infinite sigma) is the only way to reach a + // non-finite result; fall back to the old floor rather than emit -inf. + if value.is_finite() { + value + } else { + f64::MIN_POSITIVE.ln() + } } #[cfg(test)] @@ -113,16 +125,16 @@ mod tests { let mut vars = VarStore::new(); let diff = vars.alloc(Gaussian::from_ms(0.0, 6.0)); let mut f = MarginFactor::new(diff, 5.0, 1.0); - assert!(f.evidence_cached.is_none()); + assert!(f.log_evidence_cached.is_none()); f.propagate(&mut vars); - let z = f.evidence_cached.unwrap(); - // pdf(5, 0, sqrt(37)) ≈ 0.046783 - assert!((z - 0.04678300292616668).abs() < 1e-10); + let z = f.log_evidence_cached.unwrap(); + // ln pdf(5, 0, sqrt(37)) = ln(0.046783...) + assert!((z.exp() - 0.04678300292616668).abs() < 1e-10); // Subsequent propagations don't change it. f.propagate(&mut vars); - assert_eq!(f.evidence_cached.unwrap(), z); + assert_eq!(f.log_evidence_cached.unwrap(), z); } #[test] diff --git a/src/factor/trunc.rs b/src/factor/trunc.rs index 3a4184d..49ef63b 100644 --- a/src/factor/trunc.rs +++ b/src/factor/trunc.rs @@ -1,8 +1,8 @@ use crate::{ - N_INF, approx, cdf, + N_INF, approx, factor::{Factor, VarId, VarStore}, gaussian::Gaussian, - sf, + ln_interval, ln_sf, }; /// EP truncation factor on a diff variable. @@ -19,7 +19,7 @@ pub struct TruncFactor { /// Outgoing message to the diff variable (initial: `N_INF`, the EP identity). pub(crate) msg: Gaussian, /// Cached evidence (linear, not log) computed from the cavity on first propagation. - pub(crate) evidence_cached: Option, + pub(crate) log_evidence_cached: Option, } impl TruncFactor { @@ -30,7 +30,7 @@ impl TruncFactor { margin, tie, msg: N_INF, - evidence_cached: None, + log_evidence_cached: None, } } } @@ -43,8 +43,8 @@ impl TruncFactor { let marginal = vars.get(self.diff); let cavity = marginal / self.msg; - if self.evidence_cached.is_none() { - self.evidence_cached = Some(cavity_evidence(cavity, self.margin, self.tie)); + if self.log_evidence_cached.is_none() { + self.log_evidence_cached = Some(cavity_log_evidence(cavity, self.margin, self.tie)); } let trunc = approx(cavity, self.margin, self.tie); @@ -69,38 +69,34 @@ impl Factor for TruncFactor { } fn log_evidence(&self, _vars: &VarStore) -> f64 { - self.evidence_cached.unwrap_or(1.0).ln() + self.log_evidence_cached.unwrap_or(0.0) } } -/// P(diff > margin) for non-tie, P(|diff| < margin) for tie. +/// `ln P(diff > margin)` for a win, `ln P(|diff| < margin)` for a tie. /// -/// 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 { +/// Computed in log space throughout. Two earlier shapes both lost the tail: +/// `1 - cdf(..)` cancelled away every digit of an unlikely outcome, and even +/// once that was fixed the linear probability underflows to zero past about 38 +/// sigma, where clamping reported -708 nats regardless of the truth. An upset +/// is the observation a log-evidence figure exists to notice, so it has to stay +/// exact precisely where it is smallest. +fn cavity_log_evidence(diff: Gaussian, margin: f64, tie: bool) -> f64 { let (mu, sigma) = (diff.mu(), diff.sigma()); - let raw = if tie { - 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) - } + let value = if tie { + ln_interval(-margin, margin, mu, sigma) } else { - sf(margin, mu, sigma) + ln_sf(margin, mu, sigma) }; - raw.clamp(f64::MIN_POSITIVE, 1.0) + // A degenerate cavity is the only route to a non-finite result; keep the + // old floor for it rather than letting -inf poison the whole history's sum. + if value.is_finite() { + value + } else { + f64::MIN_POSITIVE.ln() + } } #[cfg(test)] @@ -131,19 +127,19 @@ mod tests { let diff = vars.alloc(Gaussian::from_ms(2.0, 3.0)); let mut f = TruncFactor::new(diff, 0.0, false); - assert!(f.evidence_cached.is_none()); + assert!(f.log_evidence_cached.is_none()); f.propagate(&mut vars); - assert!(f.evidence_cached.is_some()); - let first = f.evidence_cached.unwrap(); + assert!(f.log_evidence_cached.is_some()); + let first = f.log_evidence_cached.unwrap(); // Evidence should be P(diff > 0) for diff ~ N(2, 9) ≈ 0.748 - assert!(first > 0.7); - assert!(first < 0.8); + assert!(first.exp() > 0.7); + assert!(first.exp() < 0.8); // Subsequent propagations don't change it. f.propagate(&mut vars); - assert_eq!(f.evidence_cached.unwrap(), first); + assert_eq!(f.log_evidence_cached.unwrap(), first); } /// The defect this guards: `1 - cdf` collapsed to zero for a surprising @@ -154,7 +150,7 @@ mod tests { #[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); + let evidence = cavity_log_evidence(Gaussian::from_ms(-9.0, 1.0), 0.0, false).exp(); assert!( evidence > f64::MIN_POSITIVE, @@ -175,23 +171,48 @@ mod tests { /// 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. + /// + /// Finiteness alone is too weak a bar — the clamped version was finite too, + /// and wrong by hundreds of nats. `log_evidence_tracks_the_analytic_tail` + /// below is the assertion that actually holds this up. #[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); + let ln_e = cavity_log_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" + ln_e.is_finite() && ln_e <= 0.0, + "mu={mu} tie={tie}: log evidence {ln_e} is not a log-probability" ); } } } + /// The clamp used to floor everything past ~38 sigma at `ln(MIN_POSITIVE)` + /// = -708, however far out the real observation was. In log space the + /// answer is a polynomial and stays exact: at 1000 sigma the truth is about + /// -500_000 nats, and -708 is not a rounding error. + #[test] + fn log_evidence_tracks_the_analytic_tail() { + for mu in [-40.0f64, -60.0, -100.0, -1000.0] { + // P(diff > 0) for diff ~ N(mu, 1), mu far below zero. + let got = cavity_log_evidence(Gaussian::from_ms(mu, 1.0), 0.0, false); + + // ln Phi(mu) ~ -mu^2/2 - ln(-mu) - ln(sqrt(2 pi)) for mu << 0. + let z = -mu; + let approx = -0.5 * z * z - z.ln() - (2.0 * std::f64::consts::PI).sqrt().ln(); + + assert!( + got < f64::MIN_POSITIVE.ln(), + "mu={mu}: {got} is still stuck on the old clamp floor" + ); + assert!( + (got - approx).abs() / approx.abs() < 1e-3, + "mu={mu}: got {got}, asymptotic expectation {approx}" + ); + } + } + #[test] fn tie_evidence_uses_two_sided() { let mut vars = VarStore::new(); @@ -201,7 +222,7 @@ mod tests { f.propagate(&mut vars); // For diff ~ N(0, 4), tie=true with margin=1: P(-1 < diff < 1) ≈ 0.383 - let ev = f.evidence_cached.unwrap(); + let ev = f.log_evidence_cached.unwrap().exp(); assert!(ev > 0.35 && ev < 0.42); } diff --git a/src/game.rs b/src/game.rs index 1dbdd17..bef541b 100644 --- a/src/game.rs +++ b/src/game.rs @@ -46,8 +46,8 @@ impl DiffFactor { /// reaches. pub(crate) fn log_evidence(&self) -> f64 { match self { - Self::Trunc(f) => f.evidence_cached.unwrap_or(1.0).ln(), - Self::Margin(f) => f.evidence_cached.unwrap_or(1.0).ln(), + Self::Trunc(f) => f.log_evidence_cached.unwrap_or(0.0), + Self::Margin(f) => f.log_evidence_cached.unwrap_or(0.0), } } @@ -733,9 +733,15 @@ mod tests { let c = p[2][0]; // T1 ULP shift: mu rounds to 25.0 (was 24.999999) under natural-parameter storage. + // + // The 1e-6-place values moved when `erfc_inv`'s sign error was fixed: + // this case runs at `p_draw = 0.5`, so it goes through `compute_margin`, + // and the margin is now 8.4e-8 from the exact quantile where it was + // 1.46e-7. Verified as movement *toward* analytic truth, not a + // regression — see `erfc_inv_matches_known_quantiles`. assert_ulps_eq!(a, Gaussian::from_ms(25.0, 6.092561), epsilon = 1e-6); - assert_ulps_eq!(b, Gaussian::from_ms(33.379314, 6.483575), epsilon = 1e-6); - assert_ulps_eq!(c, Gaussian::from_ms(16.620685, 6.483575), epsilon = 1e-6); + assert_ulps_eq!(b, Gaussian::from_ms(33.379315, 6.483576), epsilon = 1e-6); + assert_ulps_eq!(c, Gaussian::from_ms(16.620685, 6.483576), epsilon = 1e-6); } #[test] diff --git a/src/lib.rs b/src/lib.rs index 40a8762..9098ed3 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -249,7 +249,15 @@ fn erfc_inv(mut y: f64) -> f64 { let t = (-2.0 * (y / 2.0).ln()).sqrt(); - let mut x = FRAC_1_SQRT_2 * ((2.30753 + t * 0.27061) / (1.0 + t * (0.99229 + t * 0.04481)) - t); + // The leading coefficient is NEGATIVE. `rational - t` is negative here, so + // a positive coefficient mirrors the starting point to `-x0` — the + // reflection of the root. Newton then has to cross the origin to get back, + // which a fixed iteration count does not manage: measured against the true + // value, `erfc_inv(0.1)` returned 1.044 instead of 1.16309, and the error + // grew as y shrank until `compute_margin` stopped being monotone in + // `p_draw` altogether. + let mut x = + -FRAC_1_SQRT_2 * ((2.30753 + t * 0.27061) / (1.0 + t * (0.99229 + t * 0.04481)) - t); for _ in 0..3 { let err = erfc(x) - y; @@ -318,6 +326,73 @@ fn erfcx(x: f64) -> f64 { } } +/// `ln` of the normal density at `x`. +/// +/// The density itself underflows to zero past about 38 sigma, and `ln` of a +/// clamped zero is -708 whatever the truth was. The log form is a polynomial: +/// it stays exact at any separation, and the values it produces (-5001 nats at +/// 100 sigma, -500001 at 1000) are perfectly representable. +pub(crate) fn ln_pdf(x: f64, mu: f64, sigma: f64) -> f64 { + let z = (x - mu) / sigma; + -(SQRT_TAU * sigma).ln() - 0.5 * z * z +} + +/// `ln P(X > x)` for `X ~ N(mu, sigma^2)`. +/// +/// In the upper tail the `exp(-z^2 / 2)` common to the tail integral is +/// factored out analytically via `erfcx`, so this never underflows — where +/// `sf(..).ln()` bottoms out at -708 once `erfc` itself reaches zero. +pub(crate) fn ln_sf(x: f64, mu: f64, sigma: f64) -> f64 { + let z = (x - mu) / sigma; + + if z > 0.0 { + // ln(0.5 * erfc(z/sqrt2)) with erfc(y) = exp(-y^2) * erfcx(y). + -std::f64::consts::LN_2 - 0.5 * z * z + erfcx(z / SQRT_2).ln() + } else { + // The mass here is at least a half; nothing to lose. + sf(x, mu, sigma).ln() + } +} + +/// `ln P(lo < X < hi)` for `X ~ N(mu, sigma^2)`. +/// +/// When the interval sits in a tail both endpoint probabilities underflow +/// together, so their difference is taken in scaled form with the shared +/// exponential factored out. When it straddles the mean nothing is small and +/// the direct difference is exact. +pub(crate) fn ln_interval(lo: f64, hi: f64, mu: f64, sigma: f64) -> f64 { + let z_lo = (lo - mu) / sigma; + let z_hi = (hi - mu) / sigma; + + if z_hi <= z_lo { + return f64::NEG_INFINITY; + } + + // Fold a lower-tail interval onto the upper tail; the normal is symmetric. + let (near, far) = if z_lo >= 0.0 { + (z_lo, z_hi) + } else if z_hi <= 0.0 { + (-z_hi, -z_lo) + } else { + // Straddles the mean: the interval holds a non-negligible share of the + // mass, so neither endpoint is near enough to 1 to cancel. + return (cdf(hi, mu, sigma) - cdf(lo, mu, sigma)) + .max(f64::MIN_POSITIVE) + .ln(); + }; + + let (a, b) = (near / SQRT_2, far / SQRT_2); + // b > a >= 0, so this ratio of exponentials is at most 1 and cannot overflow. + let scale = (a * a - b * b).exp(); + let bracket = erfcx(a) - scale * erfcx(b); + + if bracket <= 0.0 { + return f64::NEG_INFINITY; + } + + -std::f64::consts::LN_2 - a * a + bracket.ln() +} + 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(); @@ -746,6 +821,156 @@ mod tests { } } + /// `erfc_inv`'s initial guess had the wrong sign, putting Newton on the + /// mirror image of the root. Three fixed iterations could not cross back, + /// so the error grew as the argument shrank: at `p_draw = 0.99` the margin + /// came out 0.503 where the answer is 2.576. + #[test] + fn erfc_inv_matches_known_quantiles() { + // sqrt(2) * erfc_inv(1 - p) is the standard normal quantile + // Phi^-1((1 + p) / 2). + for (p, exact) in [ + (0.5f64, 0.674_489_750_196_081_7f64), + (0.9, 1.644_853_626_951_472_7), + (0.95, 1.959_963_984_540_054_2), + (0.99, 2.575_829_303_548_9), + (0.999, 3.290_526_731_491_896_4), + ] { + let got = SQRT_2 * erfc_inv(1.0 - p); + assert!( + (got - exact).abs() / exact < 1e-6, + "p={p}: got {got}, exact {exact}" + ); + } + } + + /// The draw margin must grow with the draw probability. It did not: it ran + /// 0.674 -> 1.476 -> 0.503 -> 0.982 as `p_draw` went 0.5 -> 0.9 -> 0.99 -> + /// 0.999, which is not a rounding error but a broken function. + #[test] + fn compute_margin_is_monotone_in_the_draw_probability() { + let mut previous = 0.0; + for p_draw in [ + 0.001f64, 0.01, 0.1, 0.25, 0.5, 0.75, 0.9, 0.99, 0.999, 0.9999, + ] { + let margin = compute_margin(p_draw, 1.0); + assert!( + margin > previous, + "p_draw={p_draw}: margin {margin} did not exceed {previous}" + ); + previous = margin; + } + } + + /// Round-tripping the margin back through the model's own CDF must recover + /// the draw probability it was built from. + #[test] + fn compute_margin_round_trips_through_the_cdf() { + for p_draw in [0.001f64, 0.1, 0.5, 0.9, 0.99, 0.999] { + for sd in [0.5f64, 1.0, 5.892_557] { + let margin = compute_margin(p_draw, sd); + // P(|X| < margin) for X ~ N(0, sd^2). + let recovered = 1.0 - 2.0 * cdf(-margin, 0.0, sd); + assert!( + (recovered - p_draw).abs() < 1e-6, + "p_draw={p_draw} sd={sd}: recovered {recovered}" + ); + } + } + } + + /// `ln_pdf`, `ln_sf` and `ln_interval` exist so evidence stays exact where + /// the linear forms underflow. Past ~38 sigma the linear value is zero and + /// its log is whatever floor it was clamped to. + #[test] + fn log_space_helpers_stay_exact_where_the_linear_forms_underflow() { + for z in [40.0f64, 60.0, 100.0, 1000.0] { + assert_eq!(pdf(z, 0.0, 1.0), 0.0, "pdf should underflow at {z}"); + assert_eq!(sf(z, 0.0, 1.0), 0.0, "sf should underflow at {z}"); + + let lp = ln_pdf(z, 0.0, 1.0); + let expected_lp = -(SQRT_TAU).ln() - 0.5 * z * z; + assert!( + (lp - expected_lp).abs() < 1e-9, + "ln_pdf({z}) = {lp}, expected {expected_lp}" + ); + + let ls = ln_sf(z, 0.0, 1.0); + // ln Phi(-z) ~ -z^2/2 - ln(z) - ln(sqrt(2 pi)) for large z. + let approx = -0.5 * z * z - z.ln() - SQRT_TAU.ln(); + assert!( + (ls - approx).abs() / approx.abs() < 1e-3, + "ln_sf({z}) = {ls}, asymptote {approx}" + ); + assert!( + ls < f64::MIN_POSITIVE.ln(), + "ln_sf({z}) still on the clamp floor" + ); + } + } + + /// Where nothing underflows, the log helpers must agree with the direct + /// forms exactly enough that nothing else in the crate shifts. + #[test] + fn log_space_helpers_agree_with_the_linear_forms_in_range() { + for z in [-3.0f64, -1.0, 0.0, 1.0, 2.0, 5.0, 10.0, 20.0] { + let lp = ln_pdf(z, 0.5, 2.0); + let direct_pdf = pdf(z, 0.5, 2.0); + assert!( + (lp.exp() - direct_pdf).abs() <= 1e-12 * direct_pdf, + "ln_pdf at {z}: {} vs {direct_pdf}", + lp.exp() + ); + + // Bounded by `erfc`'s ~1e-7, not tighter: above x = 2 `erfcx` uses + // a continued fraction accurate to ~1e-15, so the log path is the + // *more* accurate of the two and they part company at `erfc`'s + // error rather than at round-off. + let ls = ln_sf(z, 0.5, 2.0); + let direct = sf(z, 0.5, 2.0); + assert!( + (ls.exp() - direct).abs() <= 1e-6 * direct.max(1e-300), + "ln_sf at {z}: {} vs {direct}", + ls.exp() + ); + } + } + + #[test] + fn ln_interval_matches_the_direct_difference_when_nothing_is_small() { + for mu in [-2.0f64, 0.0, 0.5, 2.0] { + let direct = cdf(1.0, mu, 1.0) - cdf(-1.0, mu, 1.0); + let logged = ln_interval(-1.0, 1.0, mu, 1.0).exp(); + // See `log_space_helpers_agree_with_the_linear_forms_in_range`: + // the gap here is `erfc`'s own error, and the log path is the more + // accurate side of it. + assert!( + (logged - direct).abs() <= 1e-6 * direct, + "mu={mu}: {logged} vs {direct}" + ); + } + } + + /// A window far out in the tail: both endpoints underflow together, so the + /// difference has to be taken in scaled form. + #[test] + fn ln_interval_survives_a_window_deep_in_the_tail() { + for mu in [-50.0f64, -100.0, -1000.0] { + let logged = ln_interval(-1.0, 1.0, mu, 1.0); + assert!(logged.is_finite(), "mu={mu}: {logged}"); + assert!( + logged < f64::MIN_POSITIVE.ln(), + "mu={mu}: {logged} is stuck on the clamp floor" + ); + // Dominated by the near edge: ln P ~ ln Phi(-(|mu| - 1)). + let near = ln_sf(-1.0, mu, 1.0); + assert!( + (logged - near).abs() < 5.0, + "mu={mu}: {logged} strays from the near-edge tail {near}" + ); + } + } + #[test] fn test_quality() { let a = Gaussian::from_ms(25.0, 3.0);