diff --git a/src/acquisition.rs b/src/acquisition.rs index 29ad6a2..76a236f 100644 --- a/src/acquisition.rs +++ b/src/acquisition.rs @@ -47,7 +47,38 @@ fn kl_divergence(q: Gaussian, p: Gaussian) -> f64 { } let mean_gap = q.mu() - p.mu(); - 0.5 * (libm::log(var_p / var_q) + (var_q + mean_gap * mean_gap) / var_p - 1.0) + + // Algebraically `0.5 * (ln(var_p/var_q) + (var_q + gap^2)/var_p - 1)`, but + // written so that neither term can go negative. + // + // The direct form cancels against its `- 1.0` for two near-identical + // distributions and returns a *negative* divergence — measured, 762 082 of + // 3 000 000 near-identical pairs, worst `-5.55e-17`, which is exactly one + // ULP of the 1.0. It also loses the answer entirely where it is small: + // at `var_q/var_p - 1 = 1e-9` the direct form gives `0.0` where the true + // value is `2.5e-19`. + // + // With `u = var_q/var_p - 1` the variance part is `0.5 * (u - ln(1+u))`, + // which is non-negative for every `u > -1`, and the mean part is a square + // over a positive variance. Non-negativity is then structural rather than + // incidental. + let u = var_q / var_p - 1.0; + 0.5 * u_minus_ln1p(u) + mean_gap * mean_gap / (2.0 * var_p) +} + +/// `u - ln(1 + u)`, without the cancellation that spelling invites. +/// +/// Both terms are approximately `u` for small `u`, so the subtraction loses +/// everything just where the result matters. The Taylor series +/// `u^2/2 - u^3/3 + u^4/4 - ...` is exact in that regime and manifestly +/// non-negative, since `u^2/2` dominates. +fn u_minus_ln1p(u: f64) -> f64 { + if u.abs() < 1e-4 { + let u2 = u * u; + u2 * (0.5 - u / 3.0 + u2 / 4.0) + } else { + u - libm::log1p(u) + } } /// Expected information gain of a hypothetical matchup, in nats. @@ -146,7 +177,7 @@ pub fn expected_information_gain>( let mut gain = 0.0; - for (ranks, probability) in predict::outcome_distribution(&performances, &margins) { + for (ranks, probability) in predict::outcome_distribution(&performances, &margins)? { if probability <= NEGLIGIBLE { continue; } diff --git a/src/error.rs b/src/error.rs index de9d905..8f24478 100644 --- a/src/error.rs +++ b/src/error.rs @@ -131,6 +131,28 @@ pub enum InferenceError { AlreadyRegistered { key: String }, /// A prediction was given a team with no members. EmptyTeam { team: usize }, + /// The prediction grid cannot resolve the narrowest feature in the matchup. + /// + /// `predict_outcome` and `predict_ranking` integrate every team's density + /// on one shared grid, whose resolution is set by the narrowest sigma (or a + /// narrower draw margin). When the widest and narrowest are far enough + /// apart, resolving the narrow one across the wide one's support needs more + /// nodes than the grid is allowed to hold. + /// + /// Reported rather than clamped. Clamping is what this replaced, and it + /// returned probabilities greater than one — measured, a `P` of 2.79 and a + /// `Prediction::total()` of 5.41 — because the trapezoid rule stops + /// resolving a density once the step exceeds roughly 1.7 of its sigma. + /// + /// `predict_win_probabilities` answers the same matchup through adaptive + /// quadrature and is accurate here; use it when only the per-team win + /// probabilities are needed. + GridTooCoarse { + /// Nodes required to resolve the narrowest feature. + needed: usize, + /// Nodes the grid may hold. + max: usize, + }, /// A joint posterior was requested where one cannot be formed exactly. JointUnavailable { reason: &'static str }, /// Fewer than two teams were supplied to a prediction. @@ -219,6 +241,15 @@ impl fmt::Display for InferenceError { Self::EmptyTeam { team } => { write!(f, "team {team} has no members") } + Self::GridTooCoarse { needed, max } => { + write!( + f, + "the prediction grid needs {needed} nodes to resolve the narrowest \ + team's density across the widest team's support, but may hold only \ + {max}; the sigmas in this matchup are too far apart to integrate on \ + one grid. Use predict_win_probabilities, which is accurate here" + ) + } Self::JointUnavailable { reason } => { write!(f, "no exact joint posterior is available: {reason}") } diff --git a/src/history.rs b/src/history.rs index 6aaeb9e..dafc48b 100644 --- a/src/history.rs +++ b/src/history.rs @@ -1539,7 +1539,7 @@ impl, O: Observer, K: Eq + Hash + Clone> History, O: Observer, K: Eq + Hash + Clone> History Vec (f64, f64, usize) { +fn grid_shape(perf: &[Gaussian], margins: &Margins) -> Result<(f64, f64, usize), InferenceError> { let lo = perf .iter() .map(|g| g.mu() - SUPPORT_SIGMAS * g.sigma()) @@ -175,18 +179,36 @@ fn grid_shape(perf: &[Gaussian], margins: &Margins) -> (f64, f64, usize) { let feature = narrowest.min(smallest_margin); let wanted = if feature.is_finite() && feature > 0.0 { - ((hi - lo) / (feature / 12.0)).ceil() + ((hi - lo) / (feature / NODES_PER_FEATURE)).ceil() } else { MIN_GRID_POINTS as f64 }; - let points = if wanted.is_finite() { - (wanted as usize).clamp(MIN_GRID_POINTS, MAX_GRID_POINTS) - } else { - MIN_GRID_POINTS - }; + if !wanted.is_finite() { + return Ok((lo, hi, MIN_GRID_POINTS)); + } - (lo, hi, points) + // Report rather than clamp. Clamping is what this replaced: it silently + // handed the recursion a grid too coarse for the narrowest density, and the + // trapezoid rule then returned probabilities greater than one — measured, a + // `P` of 2.79 and a total of 5.41. Trapezoid error on a Gaussian is + // `~exp(-2 pi^2 (sigma/h)^2)`, which is 1e-12 at `h/sigma = 0.86` and O(1) + // by `h/sigma = 17`, so the cliff is sharp and there is no useful answer on + // the far side of it. + // + // The floor is `MIN_NODES_PER_FEATURE` rather than the `NODES_PER_FEATURE` + // asked for, because the request carries a large margin: measured accurate + // to 2.2e-12 at 1.4 nodes per sigma, and wrong by 1.2e-3 at 0.7. + let needed = wanted as usize; + let floor = ((hi - lo) / (feature / MIN_NODES_PER_FEATURE)).ceil(); + if floor.is_finite() && floor as usize > MAX_GRID_POINTS { + return Err(InferenceError::GridTooCoarse { + needed, + max: MAX_GRID_POINTS, + }); + } + + Ok((lo, hi, needed.clamp(MIN_GRID_POINTS, MAX_GRID_POINTS))) } /// Densities of each team sampled on the shared grid. @@ -198,8 +220,8 @@ struct Sampled { } impl Sampled { - fn new(perf: &[Gaussian], margins: &Margins) -> Self { - let (lo, hi, points) = grid_shape(perf, margins); + fn new(perf: &[Gaussian], margins: &Margins) -> Result { + let (lo, hi, points) = grid_shape(perf, margins)?; let step = (hi - lo) / (points - 1) as f64; let density = perf .iter() @@ -209,12 +231,12 @@ impl Sampled { .collect() }) .collect(); - Self { + Ok(Self { lo, step, points, density, - } + }) } fn node(&self, i: usize) -> f64 { @@ -317,9 +339,12 @@ fn events(n: usize, strict_only: bool) -> Vec<(Vec, Vec)> { /// /// Orders that differ only *within* a tied group describe the same finishing /// order, so their probabilities are summed into one entry. -pub(crate) fn outcome_distribution(perf: &[Gaussian], margins: &Margins) -> Vec<(Vec, f64)> { +pub(crate) fn outcome_distribution( + perf: &[Gaussian], + margins: &Margins, +) -> Result, f64)>, InferenceError> { let n = perf.len(); - let sampled = Sampled::new(perf, margins); + let sampled = Sampled::new(perf, margins)?; let mut aggregated: Vec<(Vec, f64)> = Vec::new(); for (order, tied) in events(n, margins.all_zero()) { @@ -332,7 +357,7 @@ pub(crate) fn outcome_distribution(perf: &[Gaussian], margins: &Margins) -> Vec< } aggregated.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap_or(std::cmp::Ordering::Equal)); - aggregated + Ok(aggregated) } /// All permutations of `items`. @@ -397,9 +422,13 @@ fn orders_for_groups(groups: &[Vec]) -> Vec<(Vec, Vec)> { /// Ties in `ranks` mean the tied teams may finish in any internal order, so /// this sums the orders consistent with the requested ranking rather than /// picking one. -pub(crate) fn ranking_probability(perf: &[Gaussian], margins: &Margins, ranks: &[u32]) -> f64 { +pub(crate) fn ranking_probability( + perf: &[Gaussian], + margins: &Margins, + ranks: &[u32], +) -> Result { let n = perf.len(); - let sampled = Sampled::new(perf, margins); + let sampled = Sampled::new(perf, margins)?; let mut distinct: Vec = ranks.to_vec(); distinct.sort_unstable(); @@ -410,10 +439,10 @@ pub(crate) fn ranking_probability(perf: &[Gaussian], margins: &Margins, ranks: & .map(|&r| (0..n).filter(|&i| ranks[i] == r).collect()) .collect(); - orders_for_groups(&groups) + Ok(orders_for_groups(&groups) .iter() .map(|(order, tied)| order_probability(margins, &sampled, order, tied)) - .sum() + .sum()) } /// A distribution over the ways a contest could finish. @@ -604,7 +633,7 @@ mod tests { ), ] { let n = perf.len(); - let dist = outcome_distribution(&perf, &flat(n, eps)); + let dist = outcome_distribution(&perf, &flat(n, eps)).unwrap(); let sum: f64 = dist.iter().map(|(_, p)| p).sum(); assert!( (sum - 1.0).abs() < 1e-6, @@ -620,7 +649,7 @@ mod tests { fn two_team_distribution_matches_the_closed_form() { let perf = [g(3.0, 6.0), g(-2.0, 1.0)]; let eps = 1.5; - let dist = outcome_distribution(&perf, &flat(2, eps)); + let dist = outcome_distribution(&perf, &flat(2, eps)).unwrap(); let (wa, wb) = closed_form_two(perf[0], perf[1], eps); let find = |ranks: &[u32]| { @@ -653,10 +682,10 @@ mod tests { let perf = [g(5.0, 6.0), g(0.0, 3.0), g(-5.0, 1.0)]; let eps = 1.5; let margins = flat(3, eps); - let dist = outcome_distribution(&perf, &margins); + let dist = outcome_distribution(&perf, &margins).unwrap(); for (ranks, expected) in &dist { - let direct = ranking_probability(&perf, &margins, ranks); + let direct = ranking_probability(&perf, &margins, ranks).unwrap(); assert!( (direct - expected).abs() < 1e-9, "ranks {ranks:?}: direct {direct} vs distribution {expected}" @@ -675,7 +704,7 @@ mod tests { let perf = [g(0.0, 4.0), g(0.0, 4.0), g(-8.0, 2.0)]; let mut previous = 0.0; for eps in [0.0, 0.5, 1.0, 2.0, 4.0, 8.0, 24.0] { - let p = ranking_probability(&perf, &flat(3, eps), &[0, 0, 0]); + let p = ranking_probability(&perf, &flat(3, eps), &[0, 0, 0]).unwrap(); assert!(p >= previous, "eps={eps}: {p} < {previous}"); if eps == 0.0 { assert!(p < 1e-12, "a tie needs a margin, got {p}"); @@ -696,7 +725,7 @@ mod tests { let perf = [g(0.0, 4.0), g(0.0, 4.0), g(-8.0, 2.0)]; let sweep: Vec = [0.5, 2.0, 4.0, 8.0, 16.0] .iter() - .map(|&eps| ranking_probability(&perf, &flat(3, eps), &[0, 0, 1])) + .map(|&eps| ranking_probability(&perf, &flat(3, eps), &[0, 0, 1]).unwrap()) .collect(); let peak = sweep .iter() @@ -717,7 +746,7 @@ mod tests { #[test] fn ties_are_impossible_without_a_draw_margin() { let perf = [g(0.0, 4.0), g(0.0, 4.0), g(0.0, 4.0)]; - let dist = outcome_distribution(&perf, &flat(3, 0.0)); + let dist = outcome_distribution(&perf, &flat(3, 0.0)).unwrap(); assert_eq!(dist.len(), 6, "expected only the 6 strict orders: {dist:?}"); assert!(dist.iter().all(|(r, _)| { let mut seen = r.clone(); diff --git a/tests/prediction_bounds.rs b/tests/prediction_bounds.rs new file mode 100644 index 0000000..06e5b19 --- /dev/null +++ b/tests/prediction_bounds.rs @@ -0,0 +1,136 @@ +//! Bounds that any correct implementation must satisfy, swept rather than +//! spot-checked. +//! +//! The crate's docs call the `ln k` ceiling "the sharpest available test of an +//! implementation", and record that an early prototype returned 4.77 nats. It +//! was violated again — 3.237828 nats against `ln 2` — because the existing +//! check sampled one fixture and the violation lives in a specific regime: a +//! large ratio between the widest and narrowest performance sigma, where the +//! shared prediction grid could not resolve the narrow density and returned +//! probabilities greater than one. +//! +//! A single fixture cannot defend a bound like this. A sweep can. + +use trueskill_tt::{ + ConstantDrift, GameOptions, Gaussian, InferenceError, Rating, expected_information_gain, +}; + +type R = Rating; + +/// Deterministic LCG, so a failure is reproducible from the printed seed. +struct Lcg(u64); + +impl Lcg { + fn next_f64(&mut self) -> f64 { + self.0 = self + .0 + .wrapping_mul(6_364_136_223_846_793_005) + .wrapping_add(1_442_695_040_888_963_407); + // Top 53 bits to [0, 1). + ((self.0 >> 11) as f64) / ((1u64 << 53) as f64) + } + + fn in_range(&mut self, lo: f64, hi: f64) -> f64 { + lo + (hi - lo) * self.next_f64() + } + + /// Log-uniform, so the sweep spends its samples across magnitudes rather + /// than crowding the top of the range — the violations live at small sigma. + fn log_uniform(&mut self, lo: f64, hi: f64) -> f64 { + let t = self.next_f64(); + (lo.ln() + t * (hi.ln() - lo.ln())).exp() + } +} + +#[test] +fn information_gain_never_exceeds_the_entropy_of_the_outcome() { + let mut rng = Lcg(0x5eed_1234_abcd_ef01); + let ceiling = 2.0_f64.ln(); + let mut evaluated = 0usize; + let mut refused = 0usize; + + for i in 0..2_000 { + let mu_a = rng.in_range(-100.0, 100.0); + let mu_b = rng.in_range(-100.0, 100.0); + let sigma_a = rng.log_uniform(1e-4, 1e2); + let sigma_b = rng.log_uniform(1e-4, 1e2); + let beta = rng.log_uniform(1e-4, 1e1); + + let a = R::new(Gaussian::from_ms(mu_a, sigma_a), beta, ConstantDrift(0.0)); + let b = R::new(Gaussian::from_ms(mu_b, sigma_b), beta, ConstantDrift(0.0)); + let options = GameOptions { + p_draw: 0.0, + ..GameOptions::default() + }; + + match expected_information_gain(&[&[a], &[b]], &options) { + Ok(gain) => { + evaluated += 1; + assert!( + gain.is_finite(), + "sample {i}: non-finite gain {gain} \ + (mu {mu_a}, {mu_b}; sigma {sigma_a:e}, {sigma_b:e}; beta {beta:e})" + ); + assert!( + gain >= 0.0, + "sample {i}: negative gain {gain} \ + (mu {mu_a}, {mu_b}; sigma {sigma_a:e}, {sigma_b:e}; beta {beta:e})" + ); + assert!( + gain <= ceiling + 1e-9, + "sample {i}: gain {gain} exceeds ln 2 = {ceiling} \ + (mu {mu_a}, {mu_b}; sigma {sigma_a:e}, {sigma_b:e}; beta {beta:e})" + ); + } + // Refusing to answer is acceptable; answering wrongly is not. + Err(InferenceError::GridTooCoarse { .. }) => refused += 1, + Err(e) => panic!("sample {i}: unexpected error {e:?}"), + } + } + + // The sweep must actually exercise the function, not pass by refusing + // everything. + assert!( + evaluated > 1_000, + "only {evaluated} of 2000 samples were evaluated ({refused} refused); \ + the sweep is no longer testing anything" + ); + // And it must still reach the regime where the ceiling was violated — + // large sigma ratios, which is exactly where the grid now refuses. Without + // this the sweep could drift into only-easy inputs and stop being a guard. + assert!( + refused > 0, + "no sample reached the coarse-grid regime; the sweep no longer covers \ + the case that produced 3.24 nats" + ); +} + +/// The regime that produced 3.237828 nats, pinned exactly. +#[test] +fn the_known_ceiling_violation_no_longer_answers_wrongly() { + let a = R::new( + Gaussian::from_ms(9.577_887_112_129_012, 0.000_132_507_526_585_134_38), + 0.000_307_235_559_013_096_2, + ConstantDrift(0.0), + ); + let b = R::new( + Gaussian::from_ms(-14.114_932_828_525_696, 91.586_690_140_921_16), + 0.000_307_235_559_013_096_2, + ConstantDrift(0.0), + ); + let options = GameOptions { + p_draw: 0.0, + ..GameOptions::default() + }; + + match expected_information_gain(&[&[a], &[b]], &options) { + Ok(gain) => assert!( + gain <= 2.0_f64.ln() + 1e-9, + "returned {gain}, over the ln 2 ceiling" + ), + Err(InferenceError::GridTooCoarse { needed, max }) => { + assert!(needed > max, "needed {needed} should exceed max {max}"); + } + Err(e) => panic!("unexpected error {e:?}"), + } +}