Merge api/gaussian-operators (#71)

`Gaussian` publicly implemented `Mul`, `Div`, `Add` and `Sub`. They were
the EP product, cavity and variance-space convolutions, and every one of
them lies to a reader who takes the operator at face value:

    a = N(10, 2)   b = N(4, 3)   c = N(1, 1)

    a * b        N(8.15, 1.66)   not 40
    a - b        sigma GREW, 2 -> sqrt(4 + 9)
    a * N(1, 0)  mu = NaN        "multiply by one"
    a / c        pi = -0.75      mu() prints a confident 0

The last is this crate's signature defect on a public operator. `Div` is
the cavity and can legitimately leave a negative precision, which is not
a distribution — and `mu()`/`sigma()` guard `pi <= 0` and report `0.0`
and `inf`, so it comes back as a plausible number with no panic, no
`Debug` marker and nothing to test against.

The four impls are now `pub(crate)` inherent methods that say what they
do: `ep_product`, `cavity`, `convolve`, `convolve_diff`, plus `scale`
for the one operation that genuinely is arithmetic. Nothing in a user's
workflow needed operator syntax; inference did, and it still has it.

`pi()` and `tau()` follow. Storing natural parameters is a performance
decision — it makes message passing two adds — not a contract. The
public surface is now exactly: `from_ms`, `from_mv`, `mu`, `sigma`,
`variance`, `probability_below`, `probability_above`. `from_mv` and
`variance` are promoted from `pub(crate)`; they are the honest pair for
callers who already hold a variance and should not pay a round trip
through the square root.

Four integration tests asserted bit-identity on `(pi, tau)`. They assert
it on `(mu, variance)` instead — still `assert_eq!`, still exact, and
`1/pi` and `tau/pi` are deterministic, so bit-equal natural parameters
give bit-equal moments. `a_nan_sigma_passes_through_from_ms` drops its
`|| g.pi().is_nan()` half: `sigma()` substitutes for `pi <= 0` and
`pi == inf`, so NaN survives to it only from a NaN precision.

`benches/gaussian.rs` is deleted. It timed two f64 additions through the
public operators, and keeping those public solely to feed it is the same
thing #73 objected to when a benchmark was dictating five public types.
The paths it covered are exercised by `batch` and `history_converge`
through the real call chain.

Closes #71.

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

Cargo.toml

@@ -33,10 +33,6 @@ bench = false
 name = "batch"
 harness = false
 
-[[bench]]
-name = "gaussian"
-harness = false
-
 [[bench]]
 name = "history_converge"
 harness = false

benches/gaussian.rs

@@ -1,53 +0,0 @@
-use criterion::{Criterion, criterion_group, criterion_main};
-use trueskill_tt::gaussian::Gaussian;
-
-fn benchmark_gaussian_arithmetic(criterion: &mut Criterion) {
-    // Define test Gaussians
-    let g1 = Gaussian::from_ms(25.0, 25.0 / 3.0);
-    let g2 = Gaussian::from_ms(0.0, 1.0);
-    let g3 = Gaussian::from_ms(1.0, 1.0);
-
-    // Benchmark addition
-    criterion.bench_function("Gaussian::add", |bencher| {
-        bencher.iter(|| g1 + g2);
-    });
-
-    // Benchmark subtraction
-    criterion.bench_function("Gaussian::sub", |bencher| {
-        bencher.iter(|| g1 - g3);
-    });
-
-    // Benchmark multiplication
-    criterion.bench_function("Gaussian::mul", |bencher| {
-        bencher.iter(|| g1 * g2);
-    });
-
-    // Benchmark division
-    // NOTE: numerator must have higher precision (smaller sigma) than the
-    // denominator in this representation; g2 (sigma=1) / g1 (sigma=8.33) is
-    // well-defined, whereas g1 / g2 underflows and panics in mu_sigma.
-    criterion.bench_function("Gaussian::div", |bencher| {
-        bencher.iter(|| g2 / g1);
-    });
-
-    // Benchmark natural parameter conversions
-    criterion.bench_function("Gaussian::pi", |bencher| {
-        bencher.iter(|| g1.pi());
-    });
-
-    criterion.bench_function("Gaussian::tau", |bencher| {
-        bencher.iter(|| g1.tau());
-    });
-
-    // Benchmark combined pi/tau operations (used in mul/div)
-    criterion.bench_function("Gaussian::pi_tau_combined", |bencher| {
-        bencher.iter(|| {
-            let pi = g1.pi();
-            let tau = g1.tau();
-            (pi, tau)
-        });
-    });
-}
-
-criterion_group!(benches, benchmark_gaussian_arithmetic);
-criterion_main!(benches);

src/acquisition.rs

@@ -159,7 +159,7 @@ pub fn expected_information_gain<T: Time, D: Drift<T>>(
         .iter()
         .map(|team| {
             team.iter()
-                .fold(crate::N00, |acc, rating| acc + rating.performance())
+                .fold(crate::N00, |acc, rating| acc.convolve(rating.performance()))
         })
         .collect();
 

src/factor/margin.rs

@@ -39,7 +39,7 @@ impl MarginFactor {
     /// exactly; `alpha < 1.0` writes `α·new_msg + (1−α)·old_msg`.
     pub(crate) fn propagate_with_alpha(&mut self, vars: &mut VarStore, alpha: f64) -> (f64, f64) {
         let marginal = vars.get(self.diff);
-        let cavity = marginal / self.msg;
+        let cavity = marginal.cavity(self.msg);
 
         if self.log_evidence_cached.is_none() {
             self.log_evidence_cached = Some(cavity_log_evidence(cavity, self.m_obs, self.sigma));
@@ -49,7 +49,7 @@ impl MarginFactor {
         let damped = self.msg.damp_natural(new_msg, alpha);
         let old_msg = self.msg;
         self.msg = damped;
-        vars.set(self.diff, cavity * damped);
+        vars.set(self.diff, cavity.ep_product(damped));
 
         old_msg.delta(damped)
     }

src/factor/trunc.rs

@@ -41,14 +41,14 @@ impl TruncFactor {
     /// exactly; `alpha < 1.0` writes `α·new_msg + (1−α)·old_msg`.
     pub(crate) fn propagate_with_alpha(&mut self, vars: &mut VarStore, alpha: f64) -> (f64, f64) {
         let marginal = vars.get(self.diff);
-        let cavity = marginal / self.msg;
+        let cavity = marginal.cavity(self.msg);
 
         if self.log_evidence_cached.is_none() {
             self.log_evidence_cached = Some(cavity_log_evidence(cavity, self.margin, self.tie));
         }
 
         let trunc = approx(cavity, self.margin, self.tie);
-        let new_msg = trunc / cavity;
+        let new_msg = trunc.cavity(cavity);
 
         let damped = self.msg.damp_natural(new_msg, alpha);
         let old_msg = self.msg;
@@ -57,7 +57,7 @@ impl TruncFactor {
         // marginal_new = cavity * stored_msg. With alpha = 1.0 this equals
         // `trunc` (since cavity * new_msg = trunc by construction); with
         // alpha < 1.0 it reflects the partially-applied update.
-        vars.set(self.diff, cavity * damped);
+        vars.set(self.diff, cavity.ep_product(damped));
 
         old_msg.delta(damped)
     }

src/game.rs

@@ -205,7 +205,12 @@ impl<T: Time, D: Drift<T>> Game<T, D> {
         self.likelihoods
             .iter()
             .zip(self.teams.iter())
-            .map(|(l, t)| l.iter().zip(t.iter()).map(|(&l, r)| l * r.prior).collect())
+            .map(|(l, t)| {
+                l.iter()
+                    .zip(t.iter())
+                    .map(|(&l, r)| l.ep_product(r.prior))
+                    .collect()
+            })
             .collect()
     }
 
@@ -364,7 +369,7 @@ impl<'a, T: Time, D: Drift<T>> GameRef<'a, T, D> {
                 .iter()
                 .zip(self.weights[t].iter())
                 .fold(N00, |p, (competitor, &w)| {
-                    p + (competitor.performance() * w)
+                    p.convolve(competitor.performance().scale(w))
                 })
         }));
 
@@ -384,28 +389,28 @@ impl<'a, T: Time, D: Drift<T>> GameRef<'a, T, D> {
             step = (0.0_f64, 0.0_f64);
 
             for (e, lf) in links[..n_diffs.saturating_sub(1)].iter_mut().enumerate() {
-                let pw = arena.team_prior[e] * arena.lhood_lose[e];
+                let pw = arena.team_prior[e].ep_product(arena.lhood_lose[e]);
-                let pl = arena.team_prior[e + 1] * arena.lhood_win[e + 1];
+                let pl = arena.team_prior[e + 1].ep_product(arena.lhood_win[e + 1]);
-                let raw = pw - pl;
+                let raw = pw.convolve_diff(pl);
-                arena.vars.set(lf.diff(), raw * lf.msg());
+                arena.vars.set(lf.diff(), raw.ep_product(lf.msg()));
                 let d = lf.propagate(&mut arena.vars, alpha);
                 step = tuple_max(step, d);
 
-                let new_ll = pw - lf.msg();
+                let new_ll = pw.convolve_diff(lf.msg());
                 step = tuple_max(step, arena.lhood_lose[e + 1].delta(new_ll));
                 arena.lhood_lose[e + 1] = new_ll;
             }
 
             for (rev_i, lf) in links[1..].iter_mut().rev().enumerate() {
                 let e = n_diffs - 1 - rev_i;
-                let pw = arena.team_prior[e] * arena.lhood_lose[e];
+                let pw = arena.team_prior[e].ep_product(arena.lhood_lose[e]);
-                let pl = arena.team_prior[e + 1] * arena.lhood_win[e + 1];
+                let pl = arena.team_prior[e + 1].ep_product(arena.lhood_win[e + 1]);
-                let raw = pw - pl;
+                let raw = pw.convolve_diff(pl);
-                arena.vars.set(lf.diff(), raw * lf.msg());
+                arena.vars.set(lf.diff(), raw.ep_product(lf.msg()));
                 let d = lf.propagate(&mut arena.vars, alpha);
                 step = tuple_max(step, d);
 
-                let new_lw = pl + lf.msg();
+                let new_lw = pl.convolve(lf.msg());
                 step = tuple_max(step, arena.lhood_win[e].delta(new_lw));
                 arena.lhood_win[e] = new_lw;
             }
@@ -415,18 +420,21 @@ impl<'a, T: Time, D: Drift<T>> GameRef<'a, T, D> {
 
         // Special case: exactly 1 diff (2-team game); loop body was empty.
         if n_diffs == 1 {
-            let raw = (arena.team_prior[0] * arena.lhood_lose[0])
+            let raw = arena.team_prior[0]
-                - (arena.team_prior[1] * arena.lhood_win[1]);
+                .ep_product(arena.lhood_lose[0])
-            arena.vars.set(links[0].diff(), raw * links[0].msg());
+                .convolve_diff(arena.team_prior[1].ep_product(arena.lhood_win[1]));
+            arena
+                .vars
+                .set(links[0].diff(), raw.ep_product(links[0].msg()));
             links[0].propagate(&mut arena.vars, alpha);
         }
 
         // Boundary updates: close the chain at both ends.
         if n_diffs > 0 {
-            let pl1 = arena.team_prior[1] * arena.lhood_win[1];
+            let pl1 = arena.team_prior[1].ep_product(arena.lhood_win[1]);
-            arena.lhood_win[0] = pl1 + links[0].msg();
+            arena.lhood_win[0] = pl1.convolve(links[0].msg());
-            let pw_last = arena.team_prior[n_teams - 2] * arena.lhood_lose[n_teams - 2];
+            let pw_last = arena.team_prior[n_teams - 2].ep_product(arena.lhood_lose[n_teams - 2]);
-            arena.lhood_lose[n_teams - 1] = pw_last - links[n_diffs - 1].msg();
+            arena.lhood_lose[n_teams - 1] = pw_last.convolve_diff(links[n_diffs - 1].msg());
         }
 
         let log_evidence: f64 = links.iter().map(DiffFactor::log_evidence).sum();
@@ -444,7 +452,7 @@ impl<'a, T: Time, D: Drift<T>> GameRef<'a, T, D> {
             .enumerate()
             .map(|(orig_i, (competitors, weights))| {
                 let si = arena.inv_buf[orig_i];
-                let m = arena.lhood_win[si] * arena.lhood_lose[si];
+                let m = arena.lhood_win[si].ep_product(arena.lhood_lose[si]);
                 // Already folded into `team_prior` at the top of the chain,
                 // indexed by sorted position.
                 let performance = arena.team_prior[si];
@@ -452,7 +460,8 @@ impl<'a, T: Time, D: Drift<T>> GameRef<'a, T, D> {
                     .iter()
                     .zip(weights.iter())
                     .map(|(competitor, &w)| {
-                        ((m - performance.exclude(competitor.performance() * w)) * (1.0 / w))
+                        m.convolve_diff(performance.exclude(competitor.performance().scale(w)))
+                            .scale(1.0 / w)
                             .forget(competitor.beta.powi(2))
                     })
                     .collect::<Vec<_>>()
@@ -504,7 +513,7 @@ impl<'a, T: Time, D: Drift<T>> GameRef<'a, T, D> {
             .map(|(l, t)| {
                 l.iter()
                     .zip(t.iter())
-                    .map(|(&l, p)| l * p.prior)
+                    .map(|(&l, p)| l.ep_product(p.prior))
                     .collect::<Vec<_>>()
             })
             .collect::<Vec<_>>()

src/gaussian.rs

@@ -1,5 +1,3 @@
-use std::ops;
-
 use crate::{MU, N_INF, SIGMA};
 
 /// A Gaussian distribution stored in natural parameters.
@@ -75,11 +73,12 @@ impl Gaussian {
     /// Construct from mean and *variance*, skipping the square-root round trip.
     ///
     /// `from_ms(mu, var.sqrt())` immediately squares the root away again to
-    /// recover `pi = 1/var`. Variance-combining operations (`Add`, `Sub`,
+    /// recover `pi = 1/var`. Variance-combining operations work in variance
-    /// `exclude`, `forget`) work in variance space throughout, so they go
+    /// space throughout, so they go through here instead and never take a
-    /// through here instead and never take a root.
+    /// root. Use it whenever you already hold a variance —
+    /// [`variance`](Gaussian::variance) is its inverse.
     #[inline]
-    pub(crate) fn from_mv(mu: f64, var: f64) -> Self {
+    pub fn from_mv(mu: f64, var: f64) -> Self {
         if var == f64::INFINITY {
             Self { pi: 0.0, tau: 0.0 }
         } else if var == 0.0 {
@@ -107,8 +106,7 @@ impl Gaussian {
     /// means more certain; `0.0` is an improper, uninformative message and
     /// `inf` is a point mass.
     #[inline]
-    #[must_use]
+    pub(crate) fn pi(&self) -> f64 {
-    pub fn pi(&self) -> f64 {
         self.pi
     }
 
@@ -117,8 +115,7 @@ impl Gaussian {
     /// Stored rather than derived, for the same reason as [`Gaussian::pi`].
     /// Meaningful only alongside `pi`: on its own it is not a location.
     #[inline]
-    #[must_use]
+    pub(crate) fn tau(&self) -> f64 {
-    pub fn tau(&self) -> f64 {
         self.tau
     }
 
@@ -142,12 +139,14 @@ impl Gaussian {
         }
     }
 
-    /// Variance, `1 / pi`, without the root-and-square of `sigma().powi(2)`.
+    /// Variance, without the root-and-square of `sigma().powi(2)`.
     ///
-    /// Mirrors `sigma()`'s treatment of the improper (`pi <= 0`) and point-mass
+    /// Mirrors [`sigma`](Gaussian::sigma)'s treatment of the improper
-    /// (`pi == inf`) cases.
+    /// (infinite) and point-mass (zero) cases, and is the inverse of
+    /// [`from_mv`](Gaussian::from_mv).
     #[inline]
-    pub(crate) fn variance(&self) -> f64 {
+    #[must_use]
+    pub fn variance(&self) -> f64 {
         if self.pi <= 0.0 {
             f64::INFINITY
         } else if self.pi.is_infinite() {
@@ -279,34 +278,65 @@ impl Default for Gaussian {
     }
 }
 
-impl ops::Add<Gaussian> for Gaussian {
+impl Gaussian {
-    type Output = Gaussian;
+    /// The EP factor **product**: multiply two messages about the same
-    /// Variance addition: (mu1 + mu2, sqrt(σ1² + σ2²)).
+    /// variable.
-    /// Used for combining performance and noise; rare relative to mul/div.
+    ///
-    fn add(self, rhs: Gaussian) -> Self::Output {
+    /// Two natural-parameter additions and no square root, which is why the
-        Self::from_mv(self.mu() + rhs.mu(), self.variance() + rhs.variance())
+    /// type stores `pi` and `tau` rather than `mu` and `sigma`. This is the
-    }
+    /// hot path.
-}
+    ///
-
+    /// Not arithmetic — `N(10, 2).ep_product(N(4, 3))` is `N(8.15, 1.66)`,
-impl ops::Sub<Gaussian> for Gaussian {
+    /// nowhere near 40. It used to be spelled `a * b`, on a public `Mul` impl,
-    type Output = Gaussian;
+    /// where that was a trap rather than a shorthand.
-    /// (mu1 - mu2, sqrt(σ1² + σ2²)). Same sigma combination as Add.
+    #[inline]
-    fn sub(self, rhs: Gaussian) -> Self::Output {
+    pub(crate) fn ep_product(self, rhs: Gaussian) -> Gaussian {
-        Self::from_mv(self.mu() - rhs.mu(), self.variance() + rhs.variance())
-    }
-}
-
-impl ops::Mul<Gaussian> for Gaussian {
-    type Output = Gaussian;
-    /// Factor product: nat-param add. Hot path — two f64 additions, no sqrt.
-    fn mul(self, rhs: Gaussian) -> Self::Output {
         Self::from_natural(self.pi + rhs.pi, self.tau + rhs.tau)
     }
-}
 
-impl ops::Mul<f64> for Gaussian {
+    /// The EP **cavity**: divide out a message this belief already absorbed.
-    type Output = Gaussian;
+    ///
-    fn mul(self, scalar: f64) -> Self::Output {
+    /// The inverse of [`ep_product`](Gaussian::ep_product), and two
+    /// subtractions rather than two additions.
+    ///
+    /// **May return an improper result.** Cancelling a message that carried
+    /// most of the precision leaves `pi <= 0`, which is not a distribution.
+    /// `mu()` reports `0.0` and `sigma()` reports `inf` for such a value —
+    /// both are the accessors' policy for "undefined", not answers. Measured:
+    /// `N(10, 2).cavity(N(1, 1))` has `pi = -0.75`, and its `mu()` prints a
+    /// confident `0`. That is why this is not a public operator.
+    #[inline]
+    pub(crate) fn cavity(self, rhs: Gaussian) -> Gaussian {
+        Self::from_natural(self.pi - rhs.pi, self.tau - rhs.tau)
+    }
+
+    /// Convolve two independent Gaussians: `N(mu1 + mu2, sqrt(v1 + v2))`.
+    ///
+    /// The distribution of a *sum* of independent variables, so the variances
+    /// add — the result is always wider than either input. Used to combine a
+    /// skill with performance noise. Goes through `from_mv` and takes no root.
+    #[inline]
+    pub(crate) fn convolve(self, rhs: Gaussian) -> Gaussian {
+        Self::from_mv(self.mu() + rhs.mu(), self.variance() + rhs.variance())
+    }
+
+    /// Convolve a *difference*: `N(mu1 - mu2, sqrt(v1 + v2))`.
+    ///
+    /// The means subtract and the variances still **add**, because a
+    /// difference of independent variables is no more certain than a sum. That
+    /// is the half that made the old `Sub` impl misleading: `a - b` grew the
+    /// sigma from 2 to `sqrt(4 + 9)`.
+    #[inline]
+    pub(crate) fn convolve_diff(self, rhs: Gaussian) -> Gaussian {
+        Self::from_mv(self.mu() - rhs.mu(), self.variance() + rhs.variance())
+    }
+
+    /// Scale by a constant: `mu` by `scalar`, `sigma` by `|scalar|`.
+    ///
+    /// The one operation that *is* ordinary arithmetic — it is the
+    /// distribution of `scalar * X`. Used for per-member weights.
+    #[inline]
+    pub(crate) fn scale(self, scalar: f64) -> Gaussian {
         if !scalar.is_finite() {
             return N_INF;
         }
@@ -321,14 +351,6 @@ impl ops::Mul<f64> for Gaussian {
     }
 }
 
-impl ops::Div<Gaussian> for Gaussian {
-    type Output = Gaussian;
-    /// Cavity: nat-param sub. Hot path — two f64 subtractions, no sqrt.
-    fn div(self, rhs: Gaussian) -> Self::Output {
-        Self::from_natural(self.pi - rhs.pi, self.tau - rhs.tau)
-    }
-}
-
 #[cfg(test)]
 mod tests {
     /// A message that did not change must report no change, even when it is
@@ -386,64 +408,64 @@ mod tests {
 
         // Subtracting such a message must not produce NaN (the original failure path).
         let proper = Gaussian::from_ms(9.75, 1.256);
-        let diff = proper - tiny_neg;
+        let diff = proper.convolve_diff(tiny_neg);
         assert!(diff.pi().is_finite() && !diff.pi().is_nan());
         assert!(diff.tau().is_finite() && !diff.tau().is_nan());
     }
 
     #[test]
-    fn test_add() {
+    fn convolve_adds_variances() {
         let n = Gaussian::from_ms(25.0, 25.0 / 3.0);
         let m = Gaussian::from_ms(0.0, 1.0);
-        let r = n + m;
+        let r = n.convolve(m);
         assert!((r.mu() - 25.0).abs() < 1e-12);
         assert!((r.sigma() - 8.393118874676116).abs() < 1e-10);
     }
 
     #[test]
-    fn test_sub() {
+    fn convolve_diff_subtracts_means_and_adds_variances() {
         let n = Gaussian::from_ms(25.0, 25.0 / 3.0);
         let m = Gaussian::from_ms(1.0, 1.0);
-        let r = n - m;
+        let r = n.convolve_diff(m);
         assert!((r.mu() - 24.0).abs() < 1e-12);
         assert!((r.sigma() - 8.393118874676116).abs() < 1e-10);
     }
 
     #[test]
-    fn test_mul() {
+    fn ep_product_is_not_arithmetic() {
         let n = Gaussian::from_ms(25.0, 25.0 / 3.0);
         let m = Gaussian::from_ms(0.0, 1.0);
-        let r = n * m;
+        let r = n.ep_product(m);
         assert!((r.mu() - 0.35488958990536273).abs() < 1e-10);
         assert!((r.sigma() - 0.992876838486922).abs() < 1e-10);
     }
 
     #[test]
-    fn test_div() {
+    fn cavity_undoes_a_product() {
         let n = Gaussian::from_ms(25.0, 25.0 / 3.0);
         let m = Gaussian::from_ms(0.0, 1.0);
-        let r = m / n;
+        let r = m.cavity(n);
         assert!((r.mu() - (-0.3652597402597402)).abs() < 1e-10);
         assert!((r.sigma() - 1.0072787050317253).abs() < 1e-10);
     }
 
     #[test]
     fn test_n00_is_add_identity() {
-        // N00 (sigma=0) is the additive identity for the variance-convolution Add op.
+        // N00 (sigma=0) is the identity for `convolve`.
-        // N_INF (sigma=inf) is the identity for the EP-product Mul op.
+        // N_INF (sigma=inf) is the identity for `ep_product`.
         let g = Gaussian::from_ms(3.0, 2.0);
         let n00 = Gaussian::from_ms(0.0, 0.0);
-        let r = n00 + g;
+        let r = n00.convolve(g);
         assert!((r.mu() - g.mu()).abs() < 1e-12);
         assert!((r.sigma() - g.sigma()).abs() < 1e-12);
     }
 
     #[test]
-    fn test_mul_is_factor_product() {
+    fn ep_product_adds_natural_parameters() {
-        // n * m in nat-params should be pi_n + pi_m, tau_n + tau_m
+        // `ep_product` in nat-params should be pi_n + pi_m, tau_n + tau_m
         let n = Gaussian::from_ms(2.0, 3.0);
         let m = Gaussian::from_ms(1.0, 2.0);
-        let r = n * m;
+        let r = n.ep_product(m);
         let expected_pi = n.pi() + m.pi();
         let expected_tau = n.tau() + m.tau();
         assert!((r.pi() - expected_pi).abs() < 1e-15);
@@ -451,10 +473,10 @@ mod tests {
     }
 
     #[test]
-    fn test_div_is_cavity() {
+    fn cavity_subtracts_natural_parameters() {
         let n = Gaussian::from_ms(2.0, 1.0);
         let m = Gaussian::from_ms(1.0, 2.0);
-        let r = n / m;
+        let r = n.cavity(m);
         let expected_pi = n.pi() - m.pi();
         let expected_tau = n.tau() - m.tau();
         assert!((r.pi() - expected_pi).abs() < 1e-15);

src/history.rs

@@ -1284,8 +1284,9 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
         let performances = skills
             .iter()
             .map(|team| {
-                team.iter()
+                team.iter().fold(crate::N00, |acc, s| {
-                    .fold(crate::N00, |acc, s| acc + s.forget(self.beta.powi(2)))
+                    acc.convolve(s.forget(self.beta.powi(2)))
+                })
             })
             .collect();
         let sizes = skills.iter().map(Vec::len).collect();

src/time_slice.rs

@@ -26,7 +26,9 @@ pub(crate) struct Skill {
 
 impl Skill {
     pub(crate) fn posterior(&self) -> Gaussian {
-        self.likelihood * self.backward * self.forward
+        self.likelihood
+            .ep_product(self.backward)
+            .ep_product(self.forward)
     }
 }
 
@@ -74,7 +76,7 @@ impl Item {
         if forward {
             Rating::new(skill.forward, r.beta, r.drift).with_drift_scale(r.drift_scale)
         } else {
-            Rating::new(skill.posterior() / self.likelihood, r.beta, r.drift)
+            Rating::new(skill.posterior().cavity(self.likelihood), r.beta, r.drift)
                 .with_drift_scale(r.drift_scale)
         }
     }
@@ -171,7 +173,7 @@ impl Event {
             for (i, item) in team.items.iter_mut().enumerate() {
                 let fresh = update.likelihoods[t][i];
                 let old_likelihood = skills.at(item.slot).likelihood;
-                let new_likelihood = (old_likelihood / item.likelihood) * fresh;
+                let new_likelihood = old_likelihood.cavity(item.likelihood).ep_product(fresh);
                 skills.at_mut(item.slot).likelihood = new_likelihood;
                 item.likelihood = fresh;
             }
@@ -435,8 +437,9 @@ impl<T: Time> TimeSlice<T> {
                 for (t, team) in event.teams.iter_mut().enumerate() {
                     for (i, item) in team.items.iter_mut().enumerate() {
                         let old_likelihood = self.skills.at(item.slot).likelihood;
-                        let new_likelihood =
+                        let new_likelihood = old_likelihood
-                            (old_likelihood / item.likelihood) * g.likelihoods[t][i];
+                            .cavity(item.likelihood)
+                            .ep_product(g.likelihoods[t][i]);
                         self.skills.at_mut(item.slot).likelihood = new_likelihood;
                         item.likelihood = g.likelihoods[t][i];
                     }
@@ -586,7 +589,7 @@ impl<T: Time> TimeSlice<T> {
 
     pub(crate) fn forward_prior_out(&self, competitor: &Index) -> Gaussian {
         let skill = self.skills.get(*competitor).unwrap();
-        skill.forward * skill.likelihood
+        skill.forward.ep_product(skill.likelihood)
     }
 
     pub(crate) fn backward_prior_out<D: Drift<T>>(
@@ -595,7 +598,7 @@ impl<T: Time> TimeSlice<T> {
         competitors: &CompetitorStore<T, D>,
     ) -> Gaussian {
         let skill = self.skills.get(*competitor).unwrap();
-        let n = skill.likelihood * skill.backward;
+        let n = skill.likelihood.ep_product(skill.backward);
         n.forget(
             competitors[*competitor]
                 .rating

tests/event_builder_members.rs

@@ -64,8 +64,11 @@ fn members_matches_the_typed_path_exactly() {
     for key in ["player", "layout_7"] {
         let a = typed.current_skill(&key).unwrap();
         let b = fluent.current_skill(&key).unwrap();
-        assert_eq!(a.pi(), b.pi(), "{key} pi");
+        // Exact equality, on the public moments rather than the natural
-        assert_eq!(a.tau(), b.tau(), "{key} tau");
+        // parameters: `mu` and `variance` are `tau/pi` and `1/pi`, so
+        // bit-equal natural parameters give bit-equal moments.
+        assert_eq!(a.mu(), b.mu(), "{key} mu");
+        assert_eq!(a.variance(), b.variance(), "{key} variance");
     }
 }
 

tests/joint_handle.rs

@@ -87,8 +87,8 @@ fn a_joint_answers_exactly_what_the_one_shot_call_does() {
         let terms = [(&a, 1.0), (&b, -1.0)];
         let one_shot = h.posterior_of(&terms).unwrap();
         let cached = joint.posterior_of(&terms).unwrap();
-        assert_eq!(one_shot.pi(), cached.pi(), "{a} - {b}");
+        assert_eq!(one_shot.mu(), cached.mu(), "{a} - {b}");
-        assert_eq!(one_shot.tau(), cached.tau(), "{a} - {b}");
+        assert_eq!(one_shot.variance(), cached.variance(), "{a} - {b}");
     }
 }
 
@@ -104,8 +104,8 @@ fn a_joint_agrees_at_a_pinned_time_too() {
             let cached = joint.posterior_of_at(time, &terms);
             match (one_shot, cached) {
                 (Ok(x), Ok(y)) => {
-                    assert_eq!(x.pi(), y.pi(), "t={time} {a} - {b}");
+                    assert_eq!(x.mu(), y.mu(), "t={time} {a} - {b}");
-                    assert_eq!(x.tau(), y.tau(), "t={time} {a} - {b}");
+                    assert_eq!(x.variance(), y.variance(), "t={time} {a} - {b}");
                 }
                 (Err(x), Err(y)) => assert_eq!(x, y, "t={time} {a} - {b}"),
                 (x, y) => panic!("t={time} {a} - {b}: disagreed on success: {x:?} vs {y:?}"),
@@ -262,8 +262,8 @@ fn unseen_competitors_match_the_one_shot_path() {
     let terms = [(&a, 1.0), (&z, -1.0)];
     let one_shot = h.posterior_of(&terms).unwrap();
     let cached = joint.posterior_of(&terms).unwrap();
-    assert_eq!(one_shot.pi(), cached.pi());
+    assert_eq!(one_shot.mu(), cached.mu());
-    assert_eq!(one_shot.tau(), cached.tau());
+    assert_eq!(one_shot.variance(), cached.variance());
 }
 
 /// A drift too small to represent must collapse, not corrupt the matrix.

tests/registration.rs

@@ -94,8 +94,8 @@ fn registering_matches_configuring_on_the_first_event() {
     };
 
     for ((k, a), (_, b)) in skills(&configured).into_iter().zip(skills(&registered)) {
-        assert_eq!(a.pi(), b.pi(), "{k} pi");
+        assert_eq!(a.mu(), b.mu(), "{k} mu");
-        assert_eq!(a.tau(), b.tau(), "{k} tau");
+        assert_eq!(a.variance(), b.variance(), "{k} variance");
     }
 }
 
@@ -225,8 +225,8 @@ fn registration_makes_the_fit_order_independent() {
     let forward = build(false);
     let backward = build(true);
     for ((k, a), (_, b)) in skills(&forward).into_iter().zip(skills(&backward)) {
-        assert_eq!(a.pi(), b.pi(), "{k} pi");
+        assert_eq!(a.mu(), b.mu(), "{k} mu");
-        assert_eq!(a.tau(), b.tau(), "{k} tau");
+        assert_eq!(a.variance(), b.variance(), "{k} variance");
     }
 }
 
@@ -246,8 +246,8 @@ fn rating_reads_back_what_was_stored() {
     .unwrap();
     let r = h.rating(&"layout").unwrap();
     assert_eq!(r.drift_scale(), 0.25);
-    assert_eq!(r.prior().pi(), PINNED.pi());
+    assert_eq!(r.prior().mu(), PINNED.mu());
-    assert_eq!(r.prior().tau(), PINNED.tau());
+    assert_eq!(r.prior().variance(), PINNED.variance());
 
     // A competitor created by an event reports the history defaults.
     h.record_winner(&"player", &"layout", 1).unwrap();

tests/validation.rs

@@ -248,9 +248,13 @@ mod builder_parameters {
         };
         let zero = fit(0.0);
         let positive = fit(25.0 / 6.0);
-        assert!(zero.pi().is_finite() && zero.pi() > 0.0);
+        // `variance` rather than `pi`: the natural parameters are the crate's
+        // internal representation and no longer public. It is the same
+        // quantity inverted, so a finite positive precision is a finite
+        // positive variance.
+        assert!(zero.variance().is_finite() && zero.variance() > 0.0);
         assert!(
-            (zero.pi() - positive.pi()).abs() > 1e-6,
+            (zero.variance() - positive.variance()).abs() > 1e-6,
             "zero beta must not merely be ignored: {zero:?} vs {positive:?}"
         );
     }
@@ -280,7 +284,10 @@ mod constructor_parameters {
     #[test]
     fn a_nan_sigma_passes_through_from_ms() {
         let g = Gaussian::from_ms(25.0, f64::NAN);
-        assert!(g.sigma().is_nan() || g.pi().is_nan());
+        // `sigma()` is NaN exactly when the precision is: it guards `pi <= 0`
+        // (reporting `inf`) and `pi == inf` (reporting `0.0`), so NaN survives
+        // only from a NaN precision.
+        assert!(g.sigma().is_nan());
     }
 
     #[test]