From 36eacf5f67f392a27ba4986976ed73a0612d8491 Mon Sep 17 00:00:00 2001 From: Anders Olsson Date: Tue, 8 Sep 2026 01:38:46 +0200 Subject: [PATCH] test: calibrate the marginals against the exact posterior MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Investigation for #46 and #47, before touching either. A scored history is linear-Gaussian, so its true joint posterior has a closed form and the crate can be checked against ground truth. Measured on five competitors: means marginal sd (crate / exact) tree (star) exact 1.000 loopy (robin) exact 0.502 On a tree the crate is exact in both. With cycles the means stay exact — the standard Gaussian-BP result, and the property ratings rely on — while marginal variances come out about half the true width. That is the opposite direction from what #47 reports, so whatever is happening in that consumer's model, the crate being conservative is not it. It also means #46 cannot be implemented as an added covariance accessor. The exact correlation between two nodes here is +0.857, so ignoring it overstates the width of a difference — but the too-narrow marginals partially cancel that, leaving 1.327x rather than 2.646x. Adding true correlations to these marginals without correcting them would give 0.765 against a true 1.524: overconfident, which is the direction the reporter specifically called unsafe. Pins the two real invariants (exactness on a tree, exact means with cycles) and deliberately only records the variance gap, since closing it is what #46 proposes. Also records the working rules this project has converged on: investigate before implementing, fix the root issue, and scout crates.io on measured accuracy rather than adoption. Refs #46, #47 Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_011hcFjNDmHXZF8URGLku5zZ --- CLAUDE.md | 14 ++ tests/marginal_calibration.rs | 258 ++++++++++++++++++++++++++++++++++ 2 files changed, 272 insertions(+) create mode 100644 tests/marginal_calibration.rs diff --git a/CLAUDE.md b/CLAUDE.md index d6c3f00..2c12113 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -24,6 +24,20 @@ is where several defects have hidden — a debug-only run is not evidence. - `approx` — `approx::AbsDiffEq` etc. for `Gaussian`. Most numerical goldens need it. - `rayon` — opt-in parallel within-slice sweep and per-slice query passes. +## Working rules + +- **Investigate before implementing.** Measure the actual behaviour first — + against an analytic reference where one exists. Several "obvious" fixes in + this repo turned out to be wrong in sign or unnecessary, and the measurement + is what caught them. +- **Fix the root issue, not the symptom.** A clamp that hides an underflow, or + a tolerance loosened to make a test pass, is a defect deferred. +- **Scout crates.io before hand-rolling numerics.** Check accuracy against an + independent reference rather than trusting downloads: `puruspe` has 1.4M + downloads and is 346 ULP off in the tail, where `libm` is 1. Fewer + dependencies is preferable, not mandatory — take the dependency when it is + measurably better. + ## Architecture A Rust port of [TrueSkillThroughTime.py](https://github.com/glandfried/TrueSkillThroughTime.py): diff --git a/tests/marginal_calibration.rs b/tests/marginal_calibration.rs new file mode 100644 index 0000000..229b601 --- /dev/null +++ b/tests/marginal_calibration.rs @@ -0,0 +1,258 @@ +//! Calibration of the crate's marginals against the EXACT posterior. +//! +//! A scored history is linear-Gaussian — `MarginFactor` encodes +//! `score_a - score_b ~ N(perf_a - perf_b, score_sigma^2)` — so the true joint +//! posterior has a closed form and the crate can be checked against ground +//! truth rather than against intuition. That is not possible for ranked +//! outcomes, whose truncation likelihood EP genuinely approximates. +//! +//! Two things are pinned here, and one is deliberately only recorded. +//! +//! **Pinned: on a tree the crate is exact**, means and variances both. Message +//! passing has no approximation to make when the factor graph has no cycles, so +//! any drift here would be a real defect. +//! +//! **Pinned: means are exact even with cycles.** This is the standard result +//! for Gaussian belief propagation (Weiss & Freeman 2001) and it is what makes +//! ratings trustworthy. +//! +//! **Recorded, not asserted: with cycles, marginal variances are too narrow.** +//! Measured on the round-robin fixture below, the crate reports sigma 1.430 +//! where the exact posterior is 2.851 — a ratio of 0.502. That is the known +//! behaviour of loopy Gaussian BP, not a bug in this crate, and it is left +//! unasserted because fixing it is exactly what #46 proposes. +//! +//! Why that matters for a consumer, and why #46 cannot be implemented as "add +//! a covariance accessor": the exact correlation between two nodes here is +//! +0.857, so a consumer computing `sqrt(sa^2 + sb^2)` for a difference +//! overstates its width. But the too-narrow marginals partially cancel that, +//! leaving 1.327x rather than 2.646x. Adding true correlations to these +//! marginals without also correcting them would give 0.765 against a true +//! 1.524 — *overconfident*, which is the unsafe direction. + +use smallvec::smallvec; +use trueskill_tt::{ConstantDrift, ConvergenceOptions, Event, History, Member, Outcome, Team}; + +const N: usize = 5; +const MU0: f64 = 0.0; +const SIGMA0: f64 = 6.0; +const BETA: f64 = 1.0; +const SCORE_SIGMA: f64 = 2.0; + +/// A STAR: every event touches c0, so the node-event graph is a tree and +/// Gaussian BP is exact. Any discrepancy here is not caused by loops. +fn tree_fixture() -> Vec<(usize, usize, f64)> { + vec![(0, 1, 3.0), (0, 2, 5.0), (0, 3, 4.0), (0, 4, 6.0)] +} + +/// (winner, loser, score_diff) +fn fixture() -> Vec<(usize, usize, f64)> { + vec![ + (0, 1, 3.0), + (0, 2, 5.0), + (1, 2, 2.0), + (3, 4, 1.0), + (0, 3, 4.0), + (1, 4, 2.5), + (2, 3, 0.5), + (0, 4, 6.0), + (1, 3, 1.5), + (2, 4, 3.0), + ] +} + +/// Invert a small symmetric positive-definite matrix by Gauss-Jordan. +fn inverse(mut a: Vec>) -> Vec> { + let n = a.len(); + let mut inv: Vec> = (0..n) + .map(|i| (0..n).map(|j| if i == j { 1.0 } else { 0.0 }).collect()) + .collect(); + for col in 0..n { + // partial pivot + let mut piv = col; + for r in col + 1..n { + if a[r][col].abs() > a[piv][col].abs() { + piv = r; + } + } + a.swap(col, piv); + inv.swap(col, piv); + let d = a[col][col]; + for j in 0..n { + a[col][j] /= d; + inv[col][j] /= d; + } + for r in 0..n { + if r == col { + continue; + } + let f = a[r][col]; + for j in 0..n { + a[r][j] -= f * a[col][j]; + inv[r][j] -= f * inv[col][j]; + } + } + } + inv +} + +/// The exact posterior of a linear-Gaussian model: +/// precision = prior precision + sum of a_k a_k^T / v_k. +fn exact_for(obs: &[(usize, usize, f64)]) -> (Vec, Vec>) { + let mut lambda = vec![vec![0.0; N]; N]; + let mut eta = vec![0.0; N]; + for (i, row) in lambda.iter_mut().enumerate() { + row[i] = 1.0 / (SIGMA0 * SIGMA0); + eta[i] = MU0 / (SIGMA0 * SIGMA0); + } + + // Each 1v1 observation: d ~ N(x_a - x_b, score_sigma^2 + 2 beta^2) + let v = SCORE_SIGMA * SCORE_SIGMA + 2.0 * BETA * BETA; + for &(a, b, d) in obs { + let mut vec_a = vec![0.0; N]; + vec_a[a] = 1.0; + vec_a[b] = -1.0; + for i in 0..N { + for j in 0..N { + lambda[i][j] += vec_a[i] * vec_a[j] / v; + } + eta[i] += vec_a[i] * d / v; + } + } + + let cov = inverse(lambda); + let mean: Vec = (0..N) + .map(|i| (0..N).map(|j| cov[i][j] * eta[j]).sum()) + .collect(); + (mean, cov) +} + +fn key(i: usize) -> &'static str { + ["c0", "c1", "c2", "c3", "c4"][i] +} + +/// Returns (worst mean error, worst sd ratio). +fn run(name: &str, obs: Vec<(usize, usize, f64)>) -> (f64, f64) { + println!("\n########## {name} ##########"); + let mut h: History = History::builder() + .mu(MU0) + .sigma(SIGMA0) + .beta(BETA) + .score_sigma(SCORE_SIGMA) + .drift(ConstantDrift(0.0)) + .convergence(ConvergenceOptions { + max_iter: 20_000, + epsilon: 1e-13, + alpha: 1.0, + }) + .build(); + + let events: Vec> = obs + .iter() + .copied() + .map(|(a, b, d)| Event { + time: 1, + teams: smallvec![ + Team::with_members([Member::new(key(a))]), + Team::with_members([Member::new(key(b))]), + ], + outcome: Outcome::scores([d, 0.0]), + }) + .collect(); + h.add_events(events).unwrap(); + let report = h.converge().unwrap(); + assert!( + report.converged, + "fixture must converge: {:?}", + report.final_step + ); + + let (mean, cov) = exact_for(&obs); + + println!("\n== marginals: crate vs the exact linear-Gaussian posterior =="); + println!( + "{:>4} {:>12} {:>12} {:>12} {:>12} {:>8}", + "node", "crate mu", "exact mu", "crate sd", "exact sd", "sd ratio" + ); + for i in 0..N { + let g = h.current_skill(&key(i)).unwrap(); + let exact_sd = cov[i][i].sqrt(); + println!( + "{:>4} {:>12.6} {:>12.6} {:>12.6} {:>12.6} {:>8.3}", + key(i), + g.mu(), + mean[i], + g.sigma(), + exact_sd, + g.sigma() / exact_sd + ); + } + + let mut worst_mean = 0.0f64; + let mut worst_ratio_gap = 0.0f64; + for i in 0..N { + let g = h.current_skill(&key(i)).unwrap(); + worst_mean = worst_mean.max((g.mu() - mean[i]).abs()); + worst_ratio_gap = worst_ratio_gap.max((g.sigma() / cov[i][i].sqrt() - 1.0).abs()); + } + + println!("\n== what a consumer actually computes for a DIFFERENCE =="); + println!( + "{:>8} {:>12} {:>14} {:>14} {:>12}", + "pair", "exact", "naive(exact)", "naive(crate)", "crate err" + ); + for i in 0..N { + for j in i + 1..N { + if i != 0 && j != 1 { + continue; + } + let gi = h.current_skill(&key(i)).unwrap(); + let gj = h.current_skill(&key(j)).unwrap(); + let exact_sd = (cov[i][i] + cov[j][j] - 2.0 * cov[i][j]).sqrt(); + let naive_exact = (cov[i][i] + cov[j][j]).sqrt(); + let naive_crate = (gi.sigma().powi(2) + gj.sigma().powi(2)).sqrt(); + let corr = cov[i][j] / (cov[i][i].sqrt() * cov[j][j].sqrt()); + println!( + "{:>8} {:>12.6} {:>14.6} {:>14.6} {:>11.3}x (corr {corr:.4})", + format!("{}-{}", key(i), key(j)), + exact_sd, + naive_exact, + naive_crate, + naive_crate / exact_sd + ); + } + } + + (worst_mean, worst_ratio_gap) +} + +/// With no cycles there is nothing for message passing to approximate. +#[test] +fn on_a_tree_the_marginals_are_exact() { + let (mean_err, sd_gap) = run("TREE (star: no loops, BP is exact)", tree_fixture()); + assert!( + mean_err < 1e-9, + "tree means should be exact, worst error {mean_err}" + ); + assert!( + sd_gap < 1e-9, + "tree sigmas should be exact, worst ratio gap {sd_gap}" + ); +} + +/// With cycles the means stay exact — the property ratings depend on — while +/// the variances do not. The variance gap is measured and reported rather than +/// asserted; see the module docs. +#[test] +fn with_cycles_the_means_stay_exact_but_the_variances_shrink() { + let (mean_err, sd_gap) = run("LOOPY (round robin)", fixture()); + assert!( + mean_err < 1e-9, + "loopy means must still be exact, worst error {mean_err}" + ); + assert!( + sd_gap > 0.1, + "the loopy variance gap is the premise of #46; if it has closed, that \ + issue and these docs need revisiting (worst ratio gap {sd_gap})" + ); +}