fix!: make the joint span slices, not just the latest one

`posterior_of` shipped in 0.5.0 reading a single slice. Measured against
a real Through-Time history that answers almost nothing: ustat's round
fit is 76 per-day slices whose last one holds a solo round, so 0 of 55
pair differences resolved and the single node that did was degenerate —
a one-competitor slice has no correlation to account for and returns the
marginal unchanged.

That was my mistake, and the fixture chose it. I validated against
single-slice histories, which is exactly the shape that cannot reveal
the problem. In a library whose premise is skill over time, competitors
are read at *their own* last appearance and those are different slices
by construction.

The joint is now time-expanded: one variable per appearance, linked by
the prior on a first appearance, the drift between consecutive ones, and
the within-slice event contrasts. Consecutive appearances with no drift
between them are the same variable rather than two joined by an infinite
precision, which keeps the matrix positive-definite when a competitor is
pinned with `drift_scale = 0`.

`posterior_of` now reads each competitor at their own latest appearance,
which is where `current_skill` reads them, so the two agree about which
posterior they describe. Adds `posterior_of_at(time, terms)` for a
comparison anchored to a moment, matching `learning_curve`'s reading.

Validated against a hand-written exact posterior for a two-competitor,
two-slice history — the precision matrix is spelled out in the test
rather than obtained from the crate, so it is an independent check
rather than a restatement. Also pinned: competitors last seen in
different slices now compare at all, means still agree with the
marginals, zero drift makes slice layout irrelevant, and more drift
widens a comparison across time.

BREAKING CHANGE: `posterior_of` and `expected_variance_reduction` now
consider the whole history rather than its latest slice, so results
change for any multi-slice history. `JointUnavailable` is now returned
when *any* slice holds ranked events, not just the last.

Refs #46, #47

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

src/history.rs

@@ -202,6 +202,19 @@ pub(crate) struct CompetitorConfig {
     drift_scale: Option<f64>,
 }
 
+/// The joint precision over a history's appearances, with the maps needed to
+/// address a competitor either at their latest appearance or at a given slice.
+struct TimeExpanded {
+    /// Row-major precision matrix over appearances.
+    lambda: Vec<f64>,
+    /// `(row, slice)` of each competitor's latest appearance.
+    latest: HashMap<Index, (usize, usize)>,
+    /// Row of each `(competitor, slice)` appearance.
+    at_slice: HashMap<(Index, usize), usize>,
+    /// Side length of `lambda`.
+    width: usize,
+}
+
 /// A linear functional resolved against one time slice.
 struct ResolvedTerms {
     /// Coefficients over the slice's own competitors, in its ordering.
@@ -765,19 +778,104 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
         Ok(crate::quality(&group_refs, self.beta))
     }
 
-    /// Resolve `terms` into a contrast over the slice's competitor order, the
+    /// The joint posterior precision over the whole history, time-expanded.
-    /// coefficients of any competitors the slice has never seen, and the mean.
     ///
-    /// An unseen competitor shares no event with the slice, so it is
+    /// A competitor's skill is not one variable but one per appearance, linked
-    /// independent of everything in it by construction; keeping those
+    /// by drift. That is the point of Through Time, and it is why a joint over
-    /// coefficients separate is what lets their variance be added rather than
+    /// a single slice answers almost nothing: competitors are each read at
-    /// solved for.
+    /// *their own* last appearance, and in a history with per-day or per-event
+    /// slices those are different slices. A 76-slice history whose last slice
+    /// holds one competitor can answer no pairwise question at all.
+    ///
+    /// Variables are `(competitor, appearance)`. Factors are the prior on a
+    /// first appearance, the drift between consecutive appearances, and the
+    /// within-slice event contrasts. Consecutive appearances with zero drift
+    /// variance are the same variable rather than two joined by an infinite
+    /// precision, which keeps the matrix positive-definite when a competitor is
+    /// pinned with `drift_scale = 0`.
+    ///
+    /// Returns the matrix, each competitor's row at its latest appearance, and
+    /// the row at each `(competitor, slice)` for time-addressed queries.
+    fn time_expanded_joint(&self) -> TimeExpanded {
+        let mut latest: HashMap<Index, (usize, usize)> = HashMap::new();
+        let mut at_slice: HashMap<(Index, usize), usize> = HashMap::new();
+        // Row of a competitor's previous appearance, and the drift variance
+        // separating it from the current one.
+        let mut previous: HashMap<Index, usize> = HashMap::new();
+        let mut drift_links: Vec<(usize, usize, f64)> = Vec::new();
+        let mut first_rows: Vec<(usize, Index)> = Vec::new();
+        let mut n = 0usize;
+
+        for (slice_idx, slice) in self.time_slices.iter().enumerate() {
+            for (agent, elapsed) in slice.appearances() {
+                let rating = &self.agents[agent].rating;
+                let row = match previous.get(&agent) {
+                    None => {
+                        let row = n;
+                        n += 1;
+                        first_rows.push((row, agent));
+                        row
+                    }
+                    Some(&prev) => {
+                        let drift = rating.drift_variance_for_elapsed(elapsed);
+                        if drift <= 0.0 {
+                            // No drift: the same latent skill, not two.
+                            prev
+                        } else {
+                            let row = n;
+                            n += 1;
+                            drift_links.push((prev, row, drift));
+                            row
+                        }
+                    }
+                };
+                previous.insert(agent, row);
+                latest.insert(agent, (row, slice_idx));
+                at_slice.insert((agent, slice_idx), row);
+            }
+        }
+
+        let mut lambda = vec![0.0; n * n];
+
+        for (row, agent) in first_rows {
+            lambda[row * n + row] += 1.0 / self.agents[agent].rating.prior.sigma().powi(2);
+        }
+        for (a, b, drift) in drift_links {
+            lambda[a * n + a] += 1.0 / drift;
+            lambda[b * n + b] += 1.0 / drift;
+            lambda[a * n + b] -= 1.0 / drift;
+            lambda[b * n + a] -= 1.0 / drift;
+        }
+        for (slice_idx, slice) in self.time_slices.iter().enumerate() {
+            for (contrast, noise) in slice.scored_contrasts(&self.agents) {
+                for (ia, ca) in &contrast {
+                    let ra = at_slice[&(*ia, slice_idx)];
+                    for (ib, cb) in &contrast {
+                        let rb = at_slice[&(*ib, slice_idx)];
+                        lambda[ra * n + rb] += ca * cb / noise;
+                    }
+                }
+            }
+        }
+
+        TimeExpanded {
+            lambda,
+            latest,
+            at_slice,
+            width: n,
+        }
+    }
+
+    /// Resolve `terms` into a contrast over the time-expanded rows, the
+    /// coefficients of competitors the history has never seen, and the mean.
+    ///
+    /// `row_for` picks which appearance of a competitor the caller means —
+    /// their latest, or the one at a given time.
     fn resolve_terms(
         &self,
         terms: &[(&K, f64)],
-        slice: &TimeSlice<T>,
-        row_of: &HashMap<Index, usize>,
         width: usize,
+        row_for: impl Fn(Index) -> Option<(usize, usize)>,
     ) -> Result<ResolvedTerms, InferenceError>
     where
         K: std::fmt::Debug,
@@ -790,16 +888,16 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
             let located = self
                 .keys
                 .get(*key)
-                .and_then(|index| row_of.get(&index).map(|row| (index, *row)));
+                .and_then(|index| row_for(index).map(|located| (index, located)));
 
             match located {
-                Some((index, row)) => {
+                Some((index, (row, slice_idx))) => {
                     contrast[row] += coefficient;
                     mean += coefficient
-                        * slice
+                        * self.time_slices[slice_idx]
                             .skills
                             .get(index)
-                            .expect("index came from this slice")
+                            .expect("row came from this slice")
                             .posterior()
                             .mu();
                 }
@@ -844,54 +942,131 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
     /// The mean is the same combination of the marginal means, which message
     /// passing already gets exactly right. Only the variance needs the joint.
     ///
+    /// # Which appearance each competitor is read at
+    ///
+    /// Each competitor is read at *their own* latest appearance, which is where
+    /// [`History::current_skill`] reads them too, so the two agree about which
+    /// posterior they describe. That matters in a Through-Time history: with
+    /// per-day or per-event slices, competitors are rarely all present in any
+    /// one of them. Use [`History::posterior_of_at`] to pin a time instead.
+    ///
     /// # Limitations
     ///
-    /// Currently exact only for a slice whose events are all scored, because a
+    /// Exact only for a history whose events are all scored, because a scored
-    /// scored likelihood is Gaussian and its factor can be rebuilt exactly. A
+    /// likelihood is Gaussian and its factor can be rebuilt exactly. A ranked
-    /// ranked outcome's truncation is approximated by EP, and reconstructing
+    /// outcome's truncation is approximated by EP, and reconstructing those
-    /// those factors needs the converged messages, which inference does not
+    /// factors needs the converged messages, which inference does not retain —
-    /// retain. Ranked slices return `JointUnavailable` rather than a plausible
+    /// so a history containing ranked events returns `JointUnavailable` rather
-    /// wrong number.
+    /// than a plausible wrong number.
+    ///
+    /// Cost is a dense solve over the history's *appearances*, not its
+    /// competitors: a competitor contributes one variable per slice it appears
+    /// in, minus any consecutive pair with no drift between them.
     ///
     /// # Errors
     ///
-    /// `UnknownKey` for a competitor absent from the latest slice, and
+    /// `UnknownKey` for a competitor the history has never seen, and
-    /// `JointUnavailable` if that slice contains ranked events or the system is
+    /// `JointUnavailable` for ranked events or a system that is not
-    /// not positive-definite.
+    /// positive-definite.
     pub fn posterior_of(&self, terms: &[(&K, f64)]) -> Result<Gaussian, InferenceError>
     where
         K: std::fmt::Debug,
     {
-        let slice = self
+        if self.time_slices.is_empty() {
-            .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 \
+                reason: "the history has no events",
-                         are not retained after convergence",
+            });
+        }
+        if !self.time_slices.iter().all(TimeSlice::all_scored) {
+            return Err(InferenceError::JointUnavailable {
+                reason: "the history contains ranked events, whose EP factors are \
+                         not retained after convergence",
             });
         }
 
-        let (order, lambda) = slice.joint_precision(&self.agents);
+        let TimeExpanded {
-        let mut row_of = HashMap::with_capacity(order.len());
+            lambda,
-        for (r, idx) in order.iter().enumerate() {
+            latest,
-            row_of.insert(*idx, r);
+            width,
+            ..
+        } = self.time_expanded_joint();
+
+        let ResolvedTerms {
+            contrast,
+            unseen,
+            mean,
+        } = self.resolve_terms(terms, width, |index| latest.get(&index).copied())?;
+
+        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 prior_var = self.sigma * self.sigma;
+        let variance: f64 = contrast.iter().zip(&z).map(|(c, z)| c * z).sum::<f64>()
+            + unseen.values().map(|c| c * c * prior_var).sum::<f64>();
+
+        Ok(Gaussian::from_mv(mean, variance))
+    }
+
+    /// Posterior of a linear combination, read as of `time`.
+    ///
+    /// Each competitor is taken at their latest appearance at or before `time`,
+    /// which is the same reading [`History::learning_curve`] gives. Use this
+    /// when a comparison must be anchored to a moment — "how did these two
+    /// stand at the end of last season" — rather than to wherever each
+    /// competitor was last seen.
+    ///
+    /// # Errors
+    ///
+    /// As [`History::posterior_of`], plus `UnknownKey` for a competitor with no
+    /// appearance at or before `time`.
+    pub fn posterior_of_at(&self, time: T, terms: &[(&K, f64)]) -> Result<Gaussian, InferenceError>
+    where
+        K: std::fmt::Debug,
+    {
+        if self.time_slices.is_empty() {
+            return Err(InferenceError::JointUnavailable {
+                reason: "the history has no events",
+            });
+        }
+        if !self.time_slices.iter().all(TimeSlice::all_scored) {
+            return Err(InferenceError::JointUnavailable {
+                reason: "the history contains ranked events, whose EP factors are \
+                         not retained after convergence",
+            });
+        }
+
+        let TimeExpanded {
+            lambda,
+            at_slice,
+            width,
+            ..
+        } = self.time_expanded_joint();
+
+        // Latest appearance at or before `time`, per competitor.
+        let mut as_of: HashMap<Index, (usize, usize)> = HashMap::new();
+        for (slice_idx, slice) in self.time_slices.iter().enumerate() {
+            if slice.time > time {
+                break;
+            }
+            for (agent, _) in slice.appearances() {
+                if let Some(row) = at_slice.get(&(agent, slice_idx)) {
+                    as_of.insert(agent, (*row, slice_idx));
+                }
+            }
         }
 
         let ResolvedTerms {
             contrast,
             unseen,
             mean,
-        } = self.resolve_terms(terms, slice, &row_of, order.len())?;
+        } = self.resolve_terms(terms, width, |index| as_of.get(&index).copied())?;
 
         let z =
             crate::joint::solve_spd(lambda, &contrast).ok_or(InferenceError::JointUnavailable {
-                reason: "the precision matrix is not positive-definite, which means \
+                reason: "the precision matrix is not positive-definite",
-                         a competitor has neither a proper prior nor any evidence",
             })?;
 
         let prior_var = self.sigma * self.sigma;
@@ -951,16 +1126,15 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
             });
         }
 
-        let slice = self
+        if self.time_slices.is_empty() {
-            .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 \
+                reason: "the history has no events",
-                         are not retained after convergence",
+            });
+        }
+        if !self.time_slices.iter().all(TimeSlice::all_scored) {
+            return Err(InferenceError::JointUnavailable {
+                reason: "the history contains ranked events, whose EP factors are \
+                         not retained after convergence",
             });
         }
 
@@ -983,14 +1157,15 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
             }
         }
 
-        let (order, lambda) = slice.joint_precision(&self.agents);
+        let TimeExpanded {
-        let mut row_of = HashMap::with_capacity(order.len());
+            lambda,
-        for (r, idx) in order.iter().enumerate() {
+            latest,
-            row_of.insert(*idx, r);
+            width,
-        }
+            ..
+        } = self.time_expanded_joint();
 
-        let target = self.resolve_terms(target, slice, &row_of, order.len())?;
+        let target = self.resolve_terms(target, width, |i| latest.get(&i).copied())?;
-        let matchup = self.resolve_terms(&matchup, slice, &row_of, order.len())?;
+        let matchup = self.resolve_terms(&matchup, width, |i| latest.get(&i).copied())?;
         let (target_contrast, target_unseen) = (target.contrast, target.unseen);
         let (matchup_contrast, matchup_unseen) = (matchup.contrast, matchup.unseen);
 
@@ -1003,7 +1178,7 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
         )?;
 
         let prior_var = self.sigma * self.sigma;
-        // Competitors outside the slice are independent, so they contribute
+        // Competitors outside the history are independent, so they contribute
         // only where the same key appears in both functionals.
         let cross_unseen: f64 = target_unseen
             .iter()

src/time_slice.rs

@@ -810,40 +810,26 @@ pub(crate) fn compute_elapsed<T: Time>(last: Option<&T>, current: &T) -> i64 {
 }
 
 impl<T: Time> TimeSlice<T> {
-    /// Precision matrix of the joint posterior over this slice's competitors.
+    /// This slice's scored event factors, as contrasts over 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
+    /// event's factor onto one competitor. So a joint has to be rebuilt from
-    /// from the factor structure rather than recovered from the messages.
+    /// 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
+    /// observed outcomes. The means are already exact, so only the second
-    /// 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.
+    /// Each entry is a contrast and the observation variance that sits on it.
-    /// Only scored events contribute their factors exactly; see the caller.
+    /// Ranked events contribute nothing: their truncation factors are EP
-    pub(crate) fn joint_precision<D: Drift<T>>(
+    /// approximations that inference does not retain.
+    pub(crate) fn scored_contrasts<D: Drift<T>>(
         &self,
         agents: &CompetitorStore<T, D>,
-    ) -> (Vec<Index>, Vec<f64>) {
+    ) -> Vec<(Vec<(Index, f64)>, f64)> {
-        let order: Vec<Index> = self.skills.keys().collect();
+        let mut out = Vec::new();
-        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 {
@@ -851,49 +837,48 @@ impl<T: Time> TimeSlice<T> {
             };
 
             // Teams best-first, matching the diff chain inference builds.
-            let mut order_idx: Vec<usize> = (0..event.teams.len()).collect();
+            let mut order: Vec<usize> = (0..event.teams.len()).collect();
-            order_idx.sort_by(|&a, &b| {
+            order.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) {
+            for pair in order.windows(2) {
                 let (hi, lo) = (pair[0], pair[1]);
-
+                let mut contrast: Vec<(Index, f64)> = Vec::new();
-                // 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 * agents[item.agent].rating.beta.powi(2);
-                        noise += w * w * beta * beta;
+                        contrast.push((item.agent, sign * w));
-                        *contrast.entry(row_of[&item.agent]).or_insert(0.0) += sign * w;
                     }
                 }
 
-                for (&i, &ci) in &contrast {
+                out.push((contrast, noise));
-                    for (&j, &cj) in &contrast {
-                        lambda[i * n + j] += ci * cj / noise;
-                    }
-                }
             }
         }
 
-        (order, lambda)
+        out
     }
 
-    /// True when every event here is scored, so `joint_precision` is exact.
+    /// True when every event here is scored, so the joint is exact.
     pub(crate) fn all_scored(&self) -> bool {
         self.events
             .iter()
             .all(|e| matches!(e.kind, EventKind::Scored { .. }))
     }
+
+    /// The competitors appearing in this slice, with the elapsed count since
+    /// each one's previous appearance.
+    pub(crate) fn appearances(&self) -> impl Iterator<Item = (Index, i64)> + '_ {
+        self.skills
+            .keys()
+            .map(|idx| (idx, self.skills.get(idx).expect("slice key").elapsed))
+    }
 }
 
 #[cfg(test)]

tests/time_expanded_joint.rs

@@ -0,0 +1,296 @@
+//! The joint must span slices, because Through Time reads each competitor at
+//! their own last appearance.
+//!
+//! The exact posterior of a multi-slice scored history is still Gaussian: the
+//! prior, the drift between appearances, and the scored likelihoods are all
+//! Gaussian. So it can be written out by hand and compared against, which is
+//! the check a single-slice fixture cannot make.
+
+use smallvec::smallvec;
+use trueskill_tt::{
+    ConstantDrift, ConvergenceOptions, Event, History, Member, Outcome, Team, UnknownKeys,
+};
+
+const SIGMA0: f64 = 6.0;
+const BETA: f64 = 1.0;
+const SCORE_SIGMA: f64 = 2.0;
+const GAMMA: f64 = 0.5;
+
+type H = History<i64, ConstantDrift, trueskill_tt::NullObserver, &'static str>;
+
+fn history(gamma: f64) -> H {
+    History::builder()
+        .mu(0.0)
+        .sigma(SIGMA0)
+        .beta(BETA)
+        .score_sigma(SCORE_SIGMA)
+        .drift(ConstantDrift(gamma))
+        .unknown_keys(UnknownKeys::Reject)
+        .convergence(ConvergenceOptions {
+            max_iter: 20_000,
+            epsilon: 1e-13,
+            alpha: 1.0,
+        })
+        .build()
+}
+
+fn duel(a: &'static str, b: &'static str, t: i64, sa: f64, sb: f64) -> Event<i64, &'static str> {
+    Event {
+        time: t,
+        teams: smallvec![
+            Team::with_members([Member::new(a)]),
+            Team::with_members([Member::new(b)]),
+        ],
+        outcome: Outcome::scores([sa, sb]),
+    }
+}
+
+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| f64::from(u8::from(i == j))).collect())
+        .collect();
+    for col in 0..n {
+        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
+}
+
+/// Two competitors, two slices ten units apart, one duel in each.
+///
+/// The exact precision is written out explicitly here rather than obtained
+/// from the crate, so this is an independent check rather than a restatement.
+/// Variables are `[a0, b0, a1, b1]`.
+#[test]
+fn a_two_slice_joint_matches_the_exact_posterior() {
+    let mut h = history(GAMMA);
+    h.add_events(vec![
+        duel("a", "b", 0, 5.0, 2.0),
+        duel("a", "b", 10, 4.0, 3.0),
+    ])
+    .unwrap();
+    let report = h.converge().unwrap();
+    assert!(report.converged, "{:?}", report.final_step);
+
+    let prior_prec = 1.0 / (SIGMA0 * SIGMA0);
+    let drift_prec = 1.0 / (10.0 * GAMMA * GAMMA);
+    let obs_prec = 1.0 / (SCORE_SIGMA * SCORE_SIGMA + 2.0 * BETA * BETA);
+
+    let mut lambda = vec![vec![0.0; 4]; 4];
+    // priors on the first appearances
+    lambda[0][0] += prior_prec;
+    lambda[1][1] += prior_prec;
+    // drift a0-a1 and b0-b1
+    for (p, q) in [(0usize, 2usize), (1, 3)] {
+        lambda[p][p] += drift_prec;
+        lambda[q][q] += drift_prec;
+        lambda[p][q] -= drift_prec;
+        lambda[q][p] -= drift_prec;
+    }
+    // one duel per slice: contrast (+1, -1) on that slice's variables
+    for (p, q) in [(0usize, 1usize), (2, 3)] {
+        lambda[p][p] += obs_prec;
+        lambda[q][q] += obs_prec;
+        lambda[p][q] -= obs_prec;
+        lambda[q][p] -= obs_prec;
+    }
+    let cov = inverse(lambda);
+
+    // The crate reads each competitor at their latest appearance: a1, b1.
+    let exact_gap = (cov[2][2] + cov[3][3] - 2.0 * cov[2][3]).sqrt();
+    let got = h.posterior_of(&[(&"a", 1.0), (&"b", -1.0)]).unwrap();
+    assert!(
+        (got.sigma() - exact_gap).abs() / exact_gap < 1e-9,
+        "difference: got {} exact {exact_gap}",
+        got.sigma()
+    );
+
+    let exact_single = cov[2][2].sqrt();
+    let got_single = h.posterior_of(&[(&"a", 1.0)]).unwrap();
+    assert!(
+        (got_single.sigma() - exact_single).abs() / exact_single < 1e-9,
+        "single node: got {} exact {exact_single}",
+        got_single.sigma()
+    );
+}
+
+/// The case that motivated this: competitors read at *different* slices, with
+/// the last slice holding only one of them. Under the old latest-slice joint
+/// this was `UnknownKey`.
+#[test]
+fn competitors_last_seen_in_different_slices_are_comparable() {
+    let mut h = history(GAMMA);
+    h.add_events(vec![
+        duel("a", "b", 0, 5.0, 2.0),
+        duel("a", "c", 10, 4.0, 3.0),
+        // the final slice holds one duel that does not involve b at all
+        duel("a", "c", 20, 6.0, 1.0),
+    ])
+    .unwrap();
+    let _ = h.converge().unwrap();
+
+    // b last appeared at time 0; a and c at time 20. All three must resolve.
+    for (x, y) in [("a", "b"), ("b", "c"), ("a", "c")] {
+        let g = h
+            .posterior_of(&[(&x, 1.0), (&y, -1.0)])
+            .unwrap_or_else(|e| panic!("{x} - {y} should resolve across slices: {e}"));
+        assert!(g.sigma() > 0.0 && g.sigma().is_finite());
+    }
+}
+
+/// The mean must agree with what message passing reports, which is exact even
+/// with cycles. Only the second moment needs the joint.
+#[test]
+fn means_agree_with_the_marginals() {
+    let mut h = history(GAMMA);
+    h.add_events(vec![
+        duel("a", "b", 0, 5.0, 2.0),
+        duel("b", "c", 5, 3.0, 1.0),
+        duel("a", "c", 10, 4.0, 2.0),
+    ])
+    .unwrap();
+    let _ = h.converge().unwrap();
+
+    for k in ["a", "b", "c"] {
+        let marginal = h.current_skill(&k).unwrap().mu();
+        let joint = h.posterior_of(&[(&k, 1.0)]).unwrap().mu();
+        assert!(
+            (marginal - joint).abs() < 1e-9,
+            "{k}: marginal {marginal}, joint {joint}"
+        );
+    }
+}
+
+/// With zero drift a competitor has one latent skill however many slices it
+/// appears in, so spreading the same events over time must not change the
+/// answer. This exercises the appearance-merging path.
+#[test]
+fn zero_drift_makes_slice_layout_irrelevant() {
+    let spread = {
+        let mut h = history(0.0);
+        h.add_events(vec![
+            duel("a", "b", 0, 5.0, 2.0),
+            duel("a", "b", 10, 4.0, 3.0),
+            duel("a", "b", 20, 6.0, 1.0),
+        ])
+        .unwrap();
+        let _ = h.converge().unwrap();
+        h.posterior_of(&[(&"a", 1.0), (&"b", -1.0)]).unwrap()
+    };
+    let together = {
+        let mut h = history(0.0);
+        h.add_events(vec![
+            duel("a", "b", 0, 5.0, 2.0),
+            duel("a", "b", 0, 4.0, 3.0),
+            duel("a", "b", 0, 6.0, 1.0),
+        ])
+        .unwrap();
+        let _ = h.converge().unwrap();
+        h.posterior_of(&[(&"a", 1.0), (&"b", -1.0)]).unwrap()
+    };
+
+    assert!(
+        (spread.sigma() - together.sigma()).abs() < 1e-9,
+        "zero drift: spread {} vs together {}",
+        spread.sigma(),
+        together.sigma()
+    );
+}
+
+/// More drift means less is carried forward from old evidence, so a comparison
+/// against a competitor last seen long ago must widen.
+#[test]
+fn drift_widens_a_comparison_across_time() {
+    let mut previous = 0.0;
+    for gamma in [0.0f64, 0.1, 0.5, 2.0] {
+        let mut h = history(gamma);
+        h.add_events(vec![
+            duel("a", "b", 0, 5.0, 2.0),
+            duel("a", "c", 100, 4.0, 3.0),
+        ])
+        .unwrap();
+        let _ = h.converge().unwrap();
+
+        // b was last seen at time 0; a at time 100.
+        let g = h.posterior_of(&[(&"a", 1.0), (&"b", -1.0)]).unwrap();
+        assert!(
+            g.sigma() > previous,
+            "gamma={gamma}: sigma {} did not exceed {previous}",
+            g.sigma()
+        );
+        previous = g.sigma();
+    }
+}
+
+/// `posterior_of_at` pins the reading to a moment, where `posterior_of` takes
+/// each competitor wherever they were last seen.
+#[test]
+fn posterior_of_at_reads_as_of_a_time() {
+    let mut h = history(GAMMA);
+    h.add_events(vec![
+        duel("a", "b", 0, 5.0, 2.0),
+        duel("a", "b", 10, 4.0, 3.0),
+        duel("a", "b", 20, 6.0, 1.0),
+    ])
+    .unwrap();
+    let _ = h.converge().unwrap();
+
+    let early = h.posterior_of_at(0, &[(&"a", 1.0), (&"b", -1.0)]).unwrap();
+    let late = h.posterior_of_at(20, &[(&"a", 1.0), (&"b", -1.0)]).unwrap();
+    let latest = h.posterior_of(&[(&"a", 1.0), (&"b", -1.0)]).unwrap();
+
+    // Asking as of the final slice is the same as asking for the latest.
+    assert!((late.mu() - latest.mu()).abs() < 1e-9);
+    assert!((late.sigma() - latest.sigma()).abs() < 1e-9);
+
+    // Reading at time 0 is a different quantity, and the smoothed estimate
+    // there is informed by everything that came after.
+    assert!(
+        (early.mu() - late.mu()).abs() > 1e-6,
+        "as-of-0 and as-of-20 should differ: {} vs {}",
+        early.mu(),
+        late.mu()
+    );
+
+    // A time before any event has nothing to read.
+    assert!(h.posterior_of_at(-1, &[(&"a", 1.0)]).is_err());
+}
+
+/// Times between slices resolve to the latest appearance at or before them.
+#[test]
+fn a_time_between_slices_reads_the_previous_appearance() {
+    let mut h = history(GAMMA);
+    h.add_events(vec![
+        duel("a", "b", 0, 5.0, 2.0),
+        duel("a", "b", 100, 4.0, 3.0),
+    ])
+    .unwrap();
+    let _ = h.converge().unwrap();
+
+    let at_zero = h.posterior_of_at(0, &[(&"a", 1.0)]).unwrap();
+    let between = h.posterior_of_at(50, &[(&"a", 1.0)]).unwrap();
+    assert!((at_zero.mu() - between.mu()).abs() < 1e-12);
+    assert!((at_zero.sigma() - between.sigma()).abs() < 1e-12);
+}