test: calibrate the marginals against the exact posterior

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) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011hcFjNDmHXZF8URGLku5zZ

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):

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<f64>>) -> Vec<Vec<f64>> {
+    let n = a.len();
+    let mut inv: Vec<Vec<f64>> = (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<f64>, Vec<Vec<f64>>) {
+    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<f64> = (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<i64, _, _, &'static str> = 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<Event<i64, &'static str>> = 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})"
+    );
+}