feat: add History::posterior_of for a linear combination of competitors

#46: every accessor returns a per-competitor marginal, and almost
nothing a consumer publishes is one competitor. Combining marginals
assumes independence, and competitors are correlated through every event
they share.

`posterior_of(&[(a, 1.0), (b, -1.0)])` returns the posterior of that
combination with the correlation intact. Validated against the exact
linear-Gaussian posterior on both a tree and a loopy fixture, for
differences and for single competitors: agreement to 1e-9 relative in
every case.

The investigation that preceded this is why it is not a covariance
accessor. Marginals from loopy message passing are about half the true
width, and ignoring correlation overstates a difference — the two errors
partially cancel, leaving 1.327x rather than 2.646x. Bolting true
correlations onto the existing marginals would have given 0.765 against
a true 1.524, which is overconfident: the direction the reporter
specifically called unsafe. Rebuilding the joint from the factor
structure fixes both at once, and a single-competitor query now returns
the exact marginal rather than the narrow one.

The precision matrix depends only on structure — who played whom, with
what weights and what noise — not on the observed outcomes, and the
means were already exact. So only the second moment is reconstructed.

Known limits, all deliberate and documented on the method:

- Latest slice only. A functional spanning times, such as "current
  versus career", needs the time-expanded joint and is not covered.
- Scored events only. A ranked outcome's truncation is EP-approximated
  and its converged factors are not retained after inference, so ranked
  slices return `JointUnavailable` rather than a plausible wrong number.
- Dense Cholesky, O(n^3) per query in the slice's competitor count:
  38.8us at 50, 5.66ms at 400, 49.1ms at 800. Fine for the sizes this
  serves today; caching the factorization per slice would make repeat
  queries O(n^2), and sparsity is the next step after that.

Refs #46, #47, #48

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011hcFjNDmHXZF8URGLku5zZ

src/error.rs

@@ -102,6 +102,8 @@ pub enum InferenceError {
     },
     /// A prediction was given a team with no members.
     EmptyTeam { team: 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.
     NotEnoughTeams { got: usize },
     /// The full outcome distribution was requested for too many teams.
@@ -168,6 +170,9 @@ impl fmt::Display for InferenceError {
             Self::EmptyTeam { team } => {
                 write!(f, "team {team} has no members")
             }
+            Self::JointUnavailable { reason } => {
+                write!(f, "no exact joint posterior is available: {reason}")
+            }
             Self::NotEnoughTeams { got } => {
                 write!(f, "prediction needs at least 2 teams, got {got}")
             }

src/history.rs

@@ -755,6 +755,95 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
         Ok(crate::quality(&group_refs, self.beta))
     }
 
+    /// Posterior of a linear combination of competitors' skills.
+    ///
+    /// `terms` pairs each competitor with its coefficient, so
+    /// `[(a, 1.0), (b, -1.0)]` is the difference `a - b` and
+    /// `[(score, 1.0), (layout, 1.0)]` is their sum.
+    ///
+    /// # Why this exists
+    ///
+    /// Every other accessor returns a per-competitor marginal, and combining
+    /// marginals assumes independence. Competitors are correlated through every
+    /// event they share — that coupling is the mechanism the model exists to
+    /// exploit — so `sqrt(sa^2 + sb^2)` overstates the width of a difference.
+    /// Measured against the exact posterior on a five-competitor round robin,
+    /// the correlation is +0.857 and the naive form is 2.6x too wide.
+    ///
+    /// The mean is the same combination of the marginal means, which message
+    /// passing already gets exactly right. Only the variance needs the joint.
+    ///
+    /// # Limitations
+    ///
+    /// Currently exact only for a slice whose events are all scored, because a
+    /// scored likelihood is Gaussian and its factor can be rebuilt exactly. A
+    /// ranked outcome's truncation is approximated by EP, and reconstructing
+    /// those factors needs the converged messages, which inference does not
+    /// retain. Ranked slices return `JointUnavailable` rather than a plausible
+    /// wrong number.
+    ///
+    /// # Errors
+    ///
+    /// `UnknownKey` for a competitor absent from the latest slice, and
+    /// `JointUnavailable` if that slice contains ranked events or the system is
+    /// not positive-definite.
+    pub fn posterior_of(&self, terms: &[(&K, f64)]) -> Result<Gaussian, InferenceError>
+    where
+        K: std::fmt::Debug,
+    {
+        let slice = self
+            .time_slices
+            .last()
+            .ok_or(InferenceError::JointUnavailable {
+                reason: "the history has no events",
+            })?;
+
+        if !slice.all_scored() {
+            return Err(InferenceError::JointUnavailable {
+                reason: "the latest slice contains ranked events, whose EP factors \
+                         are not retained after convergence",
+            });
+        }
+
+        let (order, lambda) = slice.joint_precision(&self.agents);
+        let mut row_of = HashMap::with_capacity(order.len());
+        for (r, idx) in order.iter().enumerate() {
+            row_of.insert(*idx, r);
+        }
+
+        let mut contrast = vec![0.0; order.len()];
+        let mut mean = 0.0;
+        for (member, (key, coefficient)) in terms.iter().enumerate() {
+            let index = self.keys.get(*key).ok_or(InferenceError::UnknownKey {
+                team: 0,
+                member,
+                key: format!("{key:?}"),
+            })?;
+            let row = *row_of.get(&index).ok_or(InferenceError::UnknownKey {
+                team: 0,
+                member,
+                key: format!("{key:?}"),
+            })?;
+            contrast[row] += coefficient;
+            mean += coefficient
+                * slice
+                    .skills
+                    .get(index)
+                    .expect("index came from this slice")
+                    .posterior()
+                    .mu();
+        }
+
+        let z =
+            crate::joint::solve_spd(lambda, &contrast).ok_or(InferenceError::JointUnavailable {
+                reason: "the precision matrix is not positive-definite, which means \
+                         a competitor has neither a proper prior nor any evidence",
+            })?;
+        let variance: f64 = contrast.iter().zip(&z).map(|(c, z)| c * z).sum();
+
+        Ok(Gaussian::from_mv(mean, variance))
+    }
+
     /// Expected information gain of running this matchup, in nats.
     ///
     /// Answers "which comparison should I run next" rather than "who will

src/joint.rs

@@ -0,0 +1,99 @@
+//! Posterior of a linear combination of competitors.
+//!
+//! Every accessor on `History` returns a per-competitor marginal, and almost
+//! nothing a consumer publishes is one competitor: "can we tell these two
+//! apart" is a difference, "what was this round worth" is a sum. Combining
+//! marginals means assuming the competitors are independent, and they are
+//! correlated through every event they share — which is the mechanism the model
+//! exists to exploit.
+//!
+//! Measured on a five-competitor round robin, the exact correlation is +0.857,
+//! so `sqrt(sa^2 + sb^2)` overstates the width of a difference by 2.6x.
+
+/// Solve `A z = b` for a symmetric positive-definite `A`, by Cholesky.
+///
+/// `a` is row-major and is consumed as scratch.
+///
+/// Returns `None` if the matrix is not positive-definite, which for a precision
+/// matrix means the model is improper — a competitor with no prior and no
+/// evidence.
+pub(crate) fn solve_spd(mut a: Vec<f64>, b: &[f64]) -> Option<Vec<f64>> {
+    let n = b.len();
+    debug_assert_eq!(a.len(), n * n);
+
+    // In-place Cholesky: A = L L^T, lower triangle.
+    for j in 0..n {
+        let mut d = a[j * n + j];
+        for k in 0..j {
+            d -= a[j * n + k] * a[j * n + k];
+        }
+        // Explicit rather than `!(d > 0.0)`: a NaN pivot must fail here too,
+        // and a negated comparison would let it through as "not positive".
+        if d.is_nan() || d <= 0.0 {
+            return None;
+        }
+        let d = d.sqrt();
+        a[j * n + j] = d;
+
+        for i in j + 1..n {
+            let mut s = a[i * n + j];
+            for k in 0..j {
+                s -= a[i * n + k] * a[j * n + k];
+            }
+            a[i * n + j] = s / d;
+        }
+    }
+
+    // Forward substitution, then back substitution.
+    let mut z = b.to_vec();
+    for i in 0..n {
+        let mut s = z[i];
+        for k in 0..i {
+            s -= a[i * n + k] * z[k];
+        }
+        z[i] = s / a[i * n + i];
+    }
+    for i in (0..n).rev() {
+        let mut s = z[i];
+        for k in i + 1..n {
+            s -= a[k * n + i] * z[k];
+        }
+        z[i] = s / a[i * n + i];
+    }
+
+    Some(z)
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn solves_a_known_system() {
+        // [[4, 1], [1, 3]] z = [1, 2]  =>  z = [1/11, 7/11]
+        let a = vec![4.0, 1.0, 1.0, 3.0];
+        let z = solve_spd(a, &[1.0, 2.0]).unwrap();
+        assert!((z[0] - 1.0 / 11.0).abs() < 1e-12, "{z:?}");
+        assert!((z[1] - 7.0 / 11.0).abs() < 1e-12, "{z:?}");
+    }
+
+    #[test]
+    fn recovers_the_inverse_diagonal() {
+        // A = [[2, -1, 0], [-1, 2, -1], [0, -1, 2]]; inverse diagonal is
+        // [0.75, 1.0, 0.75].
+        let a = vec![2.0, -1.0, 0.0, -1.0, 2.0, -1.0, 0.0, -1.0, 2.0];
+        for (i, expected) in [0.75, 1.0, 0.75].into_iter().enumerate() {
+            let mut e = vec![0.0; 3];
+            e[i] = 1.0;
+            let z = solve_spd(a.clone(), &e).unwrap();
+            assert!((z[i] - expected).abs() < 1e-12, "row {i}: {z:?}");
+        }
+    }
+
+    #[test]
+    fn rejects_a_non_positive_definite_matrix() {
+        // Singular: the second row is a multiple of the first.
+        let a = vec![1.0, 2.0, 2.0, 4.0];
+        assert!(solve_spd(a, &[1.0, 1.0]).is_none());
+    }
+}

src/lib.rs

@@ -123,6 +123,7 @@ mod game;
 pub mod gaussian;
 pub mod graph;
 mod history;
+mod joint;
 mod key_table;
 mod matrix;
 mod observer;

src/time_slice.rs

@@ -809,6 +809,93 @@ pub(crate) fn compute_elapsed<T: Time>(last: Option<&T>, current: &T) -> i64 {
     elapsed.max(0)
 }
 
+impl<T: Time> TimeSlice<T> {
+    /// Precision matrix of the joint posterior over this slice's competitors.
+    ///
+    /// Message passing produces per-competitor marginals and throws the
+    /// correlation away — `Item::likelihood` is already the projection of an
+    /// event's factor down onto one competitor. So the joint has to be rebuilt
+    /// from the factor structure rather than recovered from the messages.
+    ///
+    /// Usefully, a precision matrix depends only on *structure* — who played
+    /// whom, with what weights and what observation noise — and not on the
+    /// observed outcomes. The means are already exact (Gaussian belief
+    /// propagation gets those right even with cycles), so only the second
+    /// moment needs rebuilding.
+    ///
+    /// Returns the competitor order and the dense matrix in row-major order.
+    /// Only scored events contribute their factors exactly; see the caller.
+    pub(crate) fn joint_precision<D: Drift<T>>(
+        &self,
+        agents: &CompetitorStore<T, D>,
+    ) -> (Vec<Index>, Vec<f64>) {
+        let order: Vec<Index> = self.skills.keys().collect();
+        let n = order.len();
+        let mut row_of: HashMap<Index, usize> = HashMap::with_capacity(n);
+        for (r, idx) in order.iter().enumerate() {
+            row_of.insert(*idx, r);
+        }
+
+        let mut lambda = vec![0.0; n * n];
+
+        // Everything outside this slice enters as each competitor's forward and
+        // backward messages, which message passing treats as independent.
+        for (r, idx) in order.iter().enumerate() {
+            let skill = self.skills.get(*idx).expect("slice key has a skill");
+            lambda[r * n + r] += (skill.forward * skill.backward).pi();
+        }
+
+        for event in &self.events {
+            let EventKind::Scored { score_sigma } = event.kind else {
+                continue;
+            };
+
+            // Teams best-first, matching the diff chain inference builds.
+            let mut order_idx: Vec<usize> = (0..event.teams.len()).collect();
+            order_idx.sort_by(|&a, &b| {
+                event.teams[b]
+                    .output
+                    .partial_cmp(&event.teams[a].output)
+                    .unwrap_or(std::cmp::Ordering::Equal)
+            });
+
+            for pair in order_idx.windows(2) {
+                let (hi, lo) = (pair[0], pair[1]);
+
+                // Contrast vector, and the observation noise that sits on top
+                // of the skills: per-member performance noise plus the score
+                // noise itself.
+                let mut contrast: HashMap<usize, f64> = HashMap::new();
+                let mut noise = score_sigma * score_sigma;
+
+                for (team, sign) in [(hi, 1.0), (lo, -1.0)] {
+                    for (m, item) in event.teams[team].items.iter().enumerate() {
+                        let w = event.weights[team][m];
+                        let beta = agents[item.agent].rating.beta;
+                        noise += w * w * beta * beta;
+                        *contrast.entry(row_of[&item.agent]).or_insert(0.0) += sign * w;
+                    }
+                }
+
+                for (&i, &ci) in &contrast {
+                    for (&j, &cj) in &contrast {
+                        lambda[i * n + j] += ci * cj / noise;
+                    }
+                }
+            }
+        }
+
+        (order, lambda)
+    }
+
+    /// True when every event here is scored, so `joint_precision` is exact.
+    pub(crate) fn all_scored(&self) -> bool {
+        self.events
+            .iter()
+            .all(|e| matches!(e.kind, EventKind::Scored { .. }))
+    }
+}
+
 #[cfg(test)]
 mod tests {
     use approx::assert_ulps_eq;

tests/marginal_calibration.rs

@@ -132,8 +132,9 @@ fn key(i: usize) -> &'static str {
 }
 
 /// Returns (worst mean error, worst sd ratio).
-fn run(name: &str, obs: Vec<(usize, usize, f64)>) -> (f64, f64) {
+fn fitted(
-    println!("\n########## {name} ##########");
+    obs: &[(usize, usize, f64)],
+) -> History<i64, ConstantDrift, trueskill_tt::NullObserver, &'static str> {
     let mut h: History<i64, _, _, &'static str> = History::builder()
         .mu(MU0)
         .sigma(SIGMA0)
@@ -167,6 +168,13 @@ fn run(name: &str, obs: Vec<(usize, usize, f64)>) -> (f64, f64) {
         report.final_step
     );
 
+    h
+}
+
+/// Returns (worst mean error, worst sd ratio gap).
+fn run(name: &str, obs: Vec<(usize, usize, f64)>) -> (f64, f64) {
+    println!("\n########## {name} ##########");
+    let h = fitted(&obs);
     let (mean, cov) = exact_for(&obs);
 
     println!("\n== marginals: crate vs the exact linear-Gaussian posterior ==");
@@ -256,3 +264,104 @@ fn with_cycles_the_means_stay_exact_but_the_variances_shrink() {
          issue and these docs need revisiting (worst ratio gap {sd_gap})"
     );
 }
+
+/// The point of #46: `posterior_of` must reproduce the exact joint, including
+/// the correlation that marginals cannot express.
+#[test]
+fn posterior_of_matches_the_exact_joint() {
+    for (name, obs) in [("tree", tree_fixture()), ("loopy", fixture())] {
+        let h = fitted(&obs);
+        let (_, cov) = exact_for(&obs);
+
+        println!("\n== posterior_of vs exact ({name}) ==");
+        println!(
+            "{:>12} {:>14} {:>14} {:>10}",
+            "functional", "posterior_of", "exact", "ratio"
+        );
+
+        for (i, j) in [(0usize, 1usize), (0, 2), (1, 3), (2, 4)] {
+            let got = h
+                .posterior_of(&[(&key(i), 1.0), (&key(j), -1.0)])
+                .expect("scored slice should have a joint");
+            let exact_sd = (cov[i][i] + cov[j][j] - 2.0 * cov[i][j]).sqrt();
+            println!(
+                "{:>12} {:>14.6} {:>14.6} {:>10.4}",
+                format!("{}-{}", key(i), key(j)),
+                got.sigma(),
+                exact_sd,
+                got.sigma() / exact_sd
+            );
+            assert!(
+                (got.sigma() - exact_sd).abs() / exact_sd < 1e-9,
+                "{name} {}-{}: posterior_of gave {} where the exact joint is {exact_sd}",
+                key(i),
+                key(j),
+                got.sigma()
+            );
+        }
+
+        // A single competitor: this is where the loopy marginal was 2x narrow.
+        for (i, row) in cov.iter().enumerate() {
+            let got = h.posterior_of(&[(&key(i), 1.0)]).unwrap();
+            let exact_sd = row[i].sqrt();
+            assert!(
+                (got.sigma() - exact_sd).abs() / exact_sd < 1e-9,
+                "{name} {}: posterior_of gave {} where exact is {exact_sd}",
+                key(i),
+                got.sigma()
+            );
+        }
+        println!("  single-competitor marginals also exact");
+    }
+}
+
+/// Cost of the dense solve as the slice grows. Recorded, not asserted.
+#[test]
+#[ignore = "timing probe, run explicitly"]
+fn cost_scaling() {
+    use std::time::Instant;
+    for n in [50usize, 100, 200, 400, 800] {
+        let names: Vec<String> = (0..n).map(|i| format!("c{i}")).collect();
+        let mut h: History<i64, _, _, String> = History::builder_with_key()
+            .score_sigma(2.0)
+            .drift(ConstantDrift(0.0))
+            .convergence(ConvergenceOptions {
+                max_iter: 200,
+                epsilon: 1e-8,
+                alpha: 1.0,
+            })
+            .build();
+        let mut seed = 5u64;
+        let mut rnd = move || {
+            seed ^= seed << 13;
+            seed ^= seed >> 7;
+            seed ^= seed << 17;
+            seed
+        };
+        let events: Vec<Event<i64, String>> = (0..n * 4)
+            .map(|_| {
+                let a = (rnd() as usize) % n;
+                let mut b = (rnd() as usize) % n;
+                if b == a {
+                    b = (b + 1) % n;
+                }
+                Event {
+                    time: 1,
+                    teams: smallvec![
+                        Team::with_members([Member::new(names[a].clone())]),
+                        Team::with_members([Member::new(names[b].clone())]),
+                    ],
+                    outcome: Outcome::scores([1.0, 0.0]),
+                }
+            })
+            .collect();
+        h.add_events(events).unwrap();
+        let _ = h.converge().unwrap();
+
+        let t = Instant::now();
+        let g = h
+            .posterior_of(&[(&names[0], 1.0), (&names[1], -1.0)])
+            .unwrap();
+        println!("  n={n:>4}: {:>10.2?}  sigma {:.6}", t.elapsed(), g.sigma());
+    }
+}