feat!: N-team outcome prediction with draw mass, replacing the 2-team panic

`predict_outcome` asserted `teams.len() == 2` and returned `[p, 1 - p]`,
allocating no probability to a draw even with `p_draw > 0`. For a
draw-enabled model the numbers were simply wrong, at any team count.

It now returns `Result<Prediction, InferenceError>` and supports N teams.

Two algorithms, both deterministic:

- Who finishes first. Performances are independent Gaussians, so this
  separates into a one-dimensional integral per team rather than a
  multivariate orthant probability. Adaptive Gauss-Kronrod evaluates it
  to ~1e-15, matching the exact two-team closed form.
- A specific finishing order. The factor graph only constrains
  rank-adjacent teams, so a full order is a chain of local constraints,
  not a general orthant integral. That chain collapses into a sequential
  recursion over cumulative integrals: O(teams * grid) per order.

Fixed-node Gauss-Hermite is the obvious tool for the first and is a trap:
when a rival's sigma is small the CDF product becomes a step narrower
than the node spacing, and the nodes step over it. Measured 4.4e-4 off
the closed form on a mildly skewed matchup and 1.7e-2 on a small-sigma
one, while still returning something that looks like a probability.
Adaptive refinement is what makes that case safe, and
`win_probabilities_survive_a_rival_with_a_tiny_sigma` pins it down.

The acceptance test is an identity rather than a golden: the outcome
space is exhaustive and disjoint, so the probabilities sum to one. Any
drift is integration error and nothing else. Gauss-Hermite failed it at
4.4e-4; this holds to ~1e-9.

Also from #21: unknown keys are now reported rather than dropped, so a
team of strangers can no longer produce a confident-looking prediction.
`predict_quality` returns `Result` for the same reason.

BREAKING CHANGE: `predict_outcome` returns `Result<Prediction, _>`
instead of `Vec<f64>`; `predict_quality` returns `Result<f64, _>`.

Refs #21, #39

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

src/error.rs

@@ -40,6 +40,24 @@ pub enum InferenceError {
     },
     /// Negative precision: a Gaussian with `pi < 0` slipped into an API call.
     NegativePrecision { pi: f64 },
+    /// A prediction referenced a key the history has no skill for.
+    ///
+    /// Reported rather than skipped: dropping unknown keys turns a team of
+    /// strangers into a confident-looking probability about nobody.
+    UnknownKey { team: usize, member: usize },
+    /// A prediction was given a team with no members.
+    EmptyTeam { team: usize },
+    /// Fewer than two teams were supplied to a prediction.
+    NotEnoughTeams { got: usize },
+    /// The full outcome distribution was requested for too many teams.
+    ///
+    /// Each realisation sorts into exactly one (order, tie-pattern) event, so
+    /// the space holds `n! * 2^(n-1)` members — 1_920 at five teams, 23_040 at
+    /// six, 322_560 at seven. Past `max` this stops being something to
+    /// enumerate on a caller's behalf; ask for individual rankings with
+    /// `predict_ranking`, or for `predict_win_probabilities`, both of which
+    /// stay cheap at any team count.
+    TooManyTeams { got: usize, max: usize },
 }
 
 impl fmt::Display for InferenceError {
@@ -90,6 +108,25 @@ impl fmt::Display for InferenceError {
             Self::NegativePrecision { pi } => {
                 write!(f, "precision must be non-negative; got {pi}")
             }
+            Self::UnknownKey { team, member } => {
+                write!(
+                    f,
+                    "team {team}, member {member}: no skill recorded for this key"
+                )
+            }
+            Self::EmptyTeam { team } => {
+                write!(f, "team {team} has no members")
+            }
+            Self::NotEnoughTeams { got } => {
+                write!(f, "prediction needs at least 2 teams, got {got}")
+            }
+            Self::TooManyTeams { got, max } => {
+                write!(
+                    f,
+                    "the outcome distribution over {got} teams is too large to enumerate (limit {max}); \
+                     use predict_ranking or predict_win_probabilities instead"
+                )
+            }
         }
     }
 }

src/history.rs

@@ -9,6 +9,7 @@ use crate::{
     gaussian::Gaussian,
     key_table::KeyTable,
     observer::{NullObserver, Observer},
+    predict::Prediction,
     rating::Rating,
     sort_time,
     storage::CompetitorStore,
@@ -522,61 +523,213 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
             .sum()
     }
 
-    /// Draw-probability quality metric for the given teams (key slices).
+    /// Each team's performance Gaussian, and its member count.
     ///
-    /// Values range roughly [0, 1]; 1 == perfectly matched. Supports any
+    /// Performance is skill inflated by `beta`: the question a prediction
-    /// number of teams.
+    /// answers is "how will they do today", not "how good are they".
     ///
-    /// # Panics
+    /// # Errors
     ///
-    /// Panics if fewer than two teams are supplied, or if a team resolves to
+    /// `NotEnoughTeams`, `EmptyTeam`, or `UnknownKey`. Unknown keys are
-    /// no known competitors — keys absent from the history, or competitors
+    /// reported rather than dropped — silently skipping them would turn a team
-    /// with no recorded skill, are dropped, so a team of entirely-unknown
+    /// of strangers into a confident-looking prediction about nobody, which is
-    /// keys becomes empty. Use `lookup` to check keys first.
+    /// the failure this replaced.
-    pub fn predict_quality(&self, teams: &[&[&K]]) -> f64 {
+    fn performances(&self, teams: &[&[&K]]) -> Result<(Vec<Gaussian>, Vec<usize>), InferenceError> {
-        let groups: Vec<Vec<Gaussian>> = teams
+        if teams.len() < 2 {
-            .iter()
+            return Err(InferenceError::NotEnoughTeams { got: teams.len() });
-            .map(|team| {
-                team.iter()
-                    .filter_map(|k| self.keys.get(*k))
-                    .filter_map(|idx| {
-                        self.time_slices
-                            .iter()
-                            .rev()
-                            .find_map(|ts| ts.skills.get(idx).map(|s| s.posterior()))
-                    })
-                    .collect()
-            })
-            .collect();
-        let group_refs: Vec<&[Gaussian]> = groups.iter().map(|g| g.as_slice()).collect();
-        crate::quality(&group_refs, self.beta)
         }
 
-    /// 2-team win probability: returns `[P(team0 wins), P(team1 wins)]`.
+        let mut performances = Vec::with_capacity(teams.len());
+        let mut sizes = Vec::with_capacity(teams.len());
+
+        for (team_idx, team) in teams.iter().enumerate() {
+            if team.is_empty() {
+                return Err(InferenceError::EmptyTeam { team: team_idx });
+            }
+
+            let mut total = crate::N00;
+            for (member_idx, key) in team.iter().enumerate() {
+                let unknown = InferenceError::UnknownKey {
+                    team: team_idx,
+                    member: member_idx,
+                };
+                let index = self.keys.get(*key).ok_or(unknown.clone())?;
+                let skill = self
+                    .time_slices
+                    .iter()
+                    .rev()
+                    .find_map(|ts| ts.skills.get(index).map(|s| s.posterior()))
+                    .ok_or(unknown)?;
+                total = total + skill.forget(self.beta.powi(2));
+            }
+
+            performances.push(total);
+            sizes.push(team.len());
+        }
+
+        Ok((performances, sizes))
+    }
+
+    /// Draw margins per team pair.
     ///
-    /// Only two teams are supported.
+    /// Inference derives the margin per rank-adjacent pair from those two
+    /// teams' betas (`Game::likelihoods`), so prediction must too — a single
+    /// game-wide margin would describe a different model than the one that
+    /// will actually be fitted.
+    fn margins(&self, sizes: &[usize]) -> crate::predict::Margins {
+        let beta_sq = self.beta.powi(2);
+        let p_draw = self.p_draw;
+
+        crate::predict::Margins::new(sizes.len(), |i, j| {
+            if p_draw == 0.0 {
+                0.0
+            } else {
+                let sd = ((sizes[i] + sizes[j]) as f64 * beta_sq).sqrt();
+                crate::compute_margin(p_draw, sd)
+            }
+        })
+    }
+
+    /// Draw-probability quality metric for the given teams (key slices).
     ///
-    /// # Panics
+    /// Values range roughly `[0, 1]`; 1 == perfectly matched. Supports any
+    /// number of teams.
     ///
-    /// Panics if `teams.len() != 2`.
+    /// Note this answers "is this matchup *fair*", which is not the same as
-    pub fn predict_outcome(&self, teams: &[&[&K]]) -> Vec<f64> {
+    /// "is this matchup *informative*" — the two coincide for two evenly
-        assert_eq!(teams.len(), 2, "predict_outcome supports exactly 2 teams");
+    /// matched teams and diverge elsewhere.
-        let gather = |team: &[&K]| -> Gaussian {
+    ///
-            team.iter()
+    /// # Errors
-                .filter_map(|k| self.keys.get(*k))
+    ///
-                .filter_map(|idx| {
+    /// `NotEnoughTeams`, `EmptyTeam`, or `UnknownKey`.
+    pub fn predict_quality(&self, teams: &[&[&K]]) -> Result<f64, InferenceError> {
+        let mut groups: Vec<Vec<Gaussian>> = Vec::with_capacity(teams.len());
+
+        for (team_idx, team) in teams.iter().enumerate() {
+            if team.is_empty() {
+                return Err(InferenceError::EmptyTeam { team: team_idx });
+            }
+            let mut members = Vec::with_capacity(team.len());
+            for (member_idx, key) in team.iter().enumerate() {
+                let unknown = InferenceError::UnknownKey {
+                    team: team_idx,
+                    member: member_idx,
+                };
+                let index = self.keys.get(*key).ok_or(unknown.clone())?;
+                members.push(
                     self.time_slices
                         .iter()
                         .rev()
-                        .find_map(|ts| ts.skills.get(idx).map(|s| s.posterior()))
+                        .find_map(|ts| ts.skills.get(index).map(|s| s.posterior()))
-                })
+                        .ok_or(unknown)?,
-                .fold(crate::N00, |acc, g| acc + g.forget(self.beta.powi(2)))
+                );
-        };
+            }
-        let a = gather(teams[0]);
+            groups.push(members);
-        let b = gather(teams[1]);
+        }
-        let diff = a - b;
+
-        let p_a = 1.0 - crate::cdf(0.0, diff.mu(), diff.sigma());
+        if groups.len() < 2 {
-        vec![p_a, 1.0 - p_a]
+            return Err(InferenceError::NotEnoughTeams { got: groups.len() });
+        }
+
+        let group_refs: Vec<&[Gaussian]> = groups.iter().map(Vec::as_slice).collect();
+        Ok(crate::quality(&group_refs, self.beta))
+    }
+
+    /// `P(team i finishes strictly first)`, for every team.
+    ///
+    /// Supports any number of teams. Because performances are independent
+    /// Gaussians, this separates into a one-dimensional integral per team —
+    /// no multivariate orthant probability is involved — and is evaluated by
+    /// adaptive quadrature to within the precision of the underlying normal
+    /// CDF (~1e-8).
+    ///
+    /// With a zero `p_draw` these sum to one. With a positive `p_draw` the
+    /// shortfall is the probability that the top place is shared.
+    ///
+    /// Cheap at any team count: cost grows as the square of the team count,
+    /// not factorially. Prefer this to [`History::predict_outcome`] when you
+    /// only need to know who wins.
+    ///
+    /// # Errors
+    ///
+    /// `NotEnoughTeams`, `EmptyTeam`, or `UnknownKey`.
+    pub fn predict_win_probabilities(&self, teams: &[&[&K]]) -> Result<Vec<f64>, InferenceError> {
+        let (performances, sizes) = self.performances(teams)?;
+        Ok(crate::predict::win_probabilities(
+            &performances,
+            &self.margins(&sizes),
+        ))
+    }
+
+    /// The full distribution over finishing orders.
+    ///
+    /// Every entry is a rank vector — the shape [`crate::Outcome::ranking`]
+    /// takes, equal ranks meaning a tie — paired with its probability. The
+    /// entries are exhaustive and disjoint, so they sum to one; that identity
+    /// is the strongest available check on the numerics and is worth asserting
+    /// in tests via [`Prediction::total`].
+    ///
+    /// Accounts for `p_draw`: with a positive draw probability, tied outcomes
+    /// carry real mass rather than being silently omitted.
+    ///
+    /// # Cost
+    ///
+    /// This enumerates the outcome space, which holds `n! * 2^(n-1)` events —
+    /// 24 at three teams, 192 at four, 1_920 at five, 23_040 at six. Each
+    /// costs one `O(teams * grid)` pass, so this is milliseconds at three or
+    /// four teams and seconds at six. Above
+    /// [`MAX_TEAMS_FOR_DISTRIBUTION`](crate::MAX_PREDICTED_TEAMS) it returns
+    /// `TooManyTeams` rather than hanging. When you need one specific ordering
+    /// use [`History::predict_ranking`], and when you only need the winner use
+    /// [`History::predict_win_probabilities`]; both stay cheap at any size.
+    ///
+    /// # Errors
+    ///
+    /// `NotEnoughTeams`, `EmptyTeam`, `UnknownKey`, or `TooManyTeams`.
+    pub fn predict_outcome(&self, teams: &[&[&K]]) -> Result<Prediction, InferenceError> {
+        if teams.len() > crate::MAX_PREDICTED_TEAMS {
+            return Err(InferenceError::TooManyTeams {
+                got: teams.len(),
+                max: crate::MAX_PREDICTED_TEAMS,
+            });
+        }
+
+        let (performances, sizes) = self.performances(teams)?;
+        Ok(Prediction::new(crate::predict::outcome_distribution(
+            &performances,
+            &self.margins(&sizes),
+        )))
+    }
+
+    /// Probability of one specific finishing order.
+    ///
+    /// `ranks` follows [`crate::Outcome::ranking`]: lower is better, and equal
+    /// values mean those teams tied. Teams sharing a rank may finish in any
+    /// internal order, so this sums over those orders rather than picking one.
+    ///
+    /// Unlike [`History::predict_outcome`] this does not enumerate the outcome
+    /// space, so it stays cheap at any team count — use it when you know which
+    /// orderings you care about.
+    ///
+    /// # Errors
+    ///
+    /// `NotEnoughTeams`, `EmptyTeam`, `UnknownKey`, or `MismatchedShape` if
+    /// `ranks` does not have one entry per team.
+    pub fn predict_ranking(&self, teams: &[&[&K]], ranks: &[u32]) -> Result<f64, InferenceError> {
+        if ranks.len() != teams.len() {
+            return Err(InferenceError::MismatchedShape {
+                kind: "ranks vs teams",
+                expected: teams.len(),
+                got: ranks.len(),
+            });
+        }
+
+        let (performances, sizes) = self.performances(teams)?;
+        Ok(crate::predict::ranking_probability(
+            &performances,
+            &self.margins(&sizes),
+            ranks,
+        ))
     }
 
     /// Run the full forward+backward convergence loop and return a summary.

src/lib.rs

@@ -126,6 +126,8 @@ mod key_table;
 mod matrix;
 mod observer;
 mod outcome;
+mod predict;
+pub(crate) mod quadrature;
 mod rating;
 pub(crate) mod schedule;
 pub mod storage;
@@ -143,6 +145,7 @@ pub use key_table::KeyTable;
 use matrix::Matrix;
 pub use observer::{NullObserver, Observer};
 pub use outcome::Outcome;
+pub use predict::Prediction;
 pub use rating::Rating;
 pub use schedule::ScheduleReport;
 pub use time::{Time, Untimed};
@@ -155,6 +158,13 @@ pub const P_DRAW: f64 = 0.0;
 pub const EPSILON: f64 = 1e-6;
 pub const ITERATIONS: usize = 30;
 
+/// Largest team count `History::predict_outcome` will enumerate.
+///
+/// The outcome space holds `n! * 2^(n-1)` events, so it grows factorially:
+/// 1_920 at five teams, 23_040 at six, 322_560 at seven. Six is where
+/// enumerating on a caller's behalf stops being reasonable.
+pub const MAX_PREDICTED_TEAMS: usize = predict::MAX_TEAMS_FOR_DISTRIBUTION;
+
 const SQRT_TAU: f64 = 2.5066282746310002;
 
 pub const N01: Gaussian = Gaussian::from_ms(0.0, 1.0);

src/predict.rs

@@ -0,0 +1,723 @@
+//! Outcome prediction: who wins, and how likely is a given finishing order.
+//!
+//! Prediction runs on *performances*, not skills. A competitor's skill is
+//! inflated by their performance noise `beta` before any comparison, which is
+//! what separates "how good are they" from "how will they do today".
+//!
+//! Two questions, two algorithms:
+//!
+//! - **Who finishes first.** Because performances are independent Gaussians,
+//!   the probability that team `i` beats every other team separates into a
+//!   *one-dimensional* integral — no multivariate orthant integral is
+//!   involved. [`quadrature::integrate`] evaluates it to near machine
+//!   precision for a few hundred `cdf` calls.
+//! - **A specific finishing order.** The factor graph only ever constrains
+//!   rank-*adjacent* teams (see `Game::run_chain`), so the joint probability
+//!   of a full order is a chain of local constraints rather than a general
+//!   orthant probability. That chain collapses into a sequential recursion:
+//!   one cumulative integral per adjacent pair, `O(teams * grid)` overall.
+//!
+//! Both are deterministic. A sampler would have been easier to write and
+//! would have made every `predict_*` call return a slightly different number,
+//! which is not a property a rating library should have.
+
+use crate::{Gaussian, quadrature};
+
+/// Teams beyond this count make the outcome enumeration impractical.
+///
+/// Each realisation sorts into exactly one (permutation, tie-pattern) event,
+/// so the space has `n! * 2^(n-1)` members: 24 at 3 teams, 192 at 4, 1_920 at
+/// 5, 23_040 at 6. The jump to 322_560 at 7 is where enumerating stops being
+/// a reasonable thing to do on a caller's behalf.
+pub(crate) const MAX_TEAMS_FOR_DISTRIBUTION: usize = 6;
+
+/// Relative tolerance for the first-place integrals.
+///
+/// Tightening past this buys nothing: the underlying `cdf` is a rational
+/// approximation with fractional error ~1.2e-7, which contributes ~6e-9 to a
+/// finished probability and dominates any further quadrature refinement.
+const WIN_TOLERANCE: f64 = 1e-8;
+
+/// Nodes for the ranking grid, and the floor below which a grid is pointless.
+///
+/// The recursion converges as O(h^2). Measured against the exact two-team
+/// closed form, 2_048 nodes leave ~1.2e-6 of discretisation error while 8_192
+/// reach ~1e-7 — at which point the residual is the `cdf` rational
+/// approximation (~2.4e-8), not the grid, and refining further buys nothing.
+const MIN_GRID_POINTS: usize = 8_192;
+const MAX_GRID_POINTS: usize = 262_144;
+
+/// How many standard deviations of support the grid and integrals cover.
+///
+/// The normal density is below 1e-18 of its peak past nine sigma, far under
+/// the precision of everything else here.
+const SUPPORT_SIGMAS: f64 = 9.0;
+
+/// Standard normal CDF at `z`.
+fn phi(z: f64) -> f64 {
+    crate::cdf(z, 0.0, 1.0)
+}
+
+/// Normal density of `x` under `g`.
+fn density(g: Gaussian, x: f64) -> f64 {
+    let sigma = g.sigma();
+    let z = (x - g.mu()) / sigma;
+    (-0.5 * z * z).exp() / (sigma * (2.0 * std::f64::consts::PI).sqrt())
+}
+
+/// Per-pair draw margins.
+///
+/// The margin is *not* a single number for the whole game: inference derives
+/// it per rank-adjacent pair from those two teams' betas (`Game::likelihoods`).
+/// Prediction has to use the same per-pair values or it answers a question
+/// about a different model than the one that will actually be fitted.
+pub(crate) struct Margins {
+    n: usize,
+    values: Vec<f64>,
+}
+
+impl Margins {
+    /// Build from a per-pair margin function.
+    pub(crate) fn new<F: Fn(usize, usize) -> f64>(n: usize, f: F) -> Self {
+        let mut values = vec![0.0; n * n];
+        for i in 0..n {
+            for j in 0..n {
+                if i != j {
+                    values[i * n + j] = f(i, j);
+                }
+            }
+        }
+        Self { n, values }
+    }
+
+    fn get(&self, i: usize, j: usize) -> f64 {
+        self.values[i * self.n + j]
+    }
+
+    /// True when no pair can draw, so every tie has probability zero.
+    fn all_zero(&self) -> bool {
+        self.values.iter().all(|&v| v == 0.0)
+    }
+}
+
+/// `P(team i finishes strictly first)` for every team.
+///
+/// Strictly means beating each rival by more than that pair's draw margin, so
+/// with a non-zero margin these sum to less than one; the shortfall is the
+/// probability that the top place is shared.
+pub(crate) fn win_probabilities(perf: &[Gaussian], margins: &Margins) -> Vec<f64> {
+    (0..perf.len())
+        .map(|i| {
+            let (mu, sigma) = (perf[i].mu(), perf[i].sigma());
+            let (lo, hi) = (mu - SUPPORT_SIGMAS * sigma, mu + SUPPORT_SIGMAS * sigma);
+
+            // Each rival's CDF turns over near its own mean plus the margin.
+            // Seeding there is what keeps a rival with a tiny sigma — a step
+            // function in disguise — from being stepped over.
+            let mut seeds = Vec::with_capacity(3 * perf.len());
+            for (j, rival) in perf.iter().enumerate().filter(|&(j, _)| j != i) {
+                let centre = rival.mu() + margins.get(i, j);
+                seeds.extend_from_slice(&[centre - rival.sigma(), centre, centre + rival.sigma()]);
+            }
+
+            quadrature::integrate(
+                |x| {
+                    let d = density(perf[i], x);
+                    if d == 0.0 {
+                        return 0.0;
+                    }
+                    let beaten: f64 = (0..perf.len())
+                        .filter(|&j| j != i)
+                        .map(|j| phi((x - margins.get(i, j) - perf[j].mu()) / perf[j].sigma()))
+                        .product();
+                    d * beaten
+                },
+                lo,
+                hi,
+                &seeds,
+                WIN_TOLERANCE,
+            )
+        })
+        .collect()
+}
+
+/// Grid bounds and resolution covering every team's support.
+///
+/// Resolution is set by the *smallest* feature in play — the narrowest sigma,
+/// or a draw margin narrower still — because that is what the recursion has to
+/// resolve. A grid sized off the widest team would step over the narrow one.
+fn grid_shape(perf: &[Gaussian], margins: &Margins) -> (f64, f64, usize) {
+    let lo = perf
+        .iter()
+        .map(|g| g.mu() - SUPPORT_SIGMAS * g.sigma())
+        .fold(f64::INFINITY, f64::min);
+    let hi = perf
+        .iter()
+        .map(|g| g.mu() + SUPPORT_SIGMAS * g.sigma())
+        .fold(f64::NEG_INFINITY, f64::max);
+
+    let narrowest = perf
+        .iter()
+        .map(Gaussian::sigma)
+        .fold(f64::INFINITY, f64::min);
+    let smallest_margin = margins
+        .values
+        .iter()
+        .copied()
+        .filter(|&m| m > 0.0)
+        .fold(f64::INFINITY, f64::min);
+
+    let feature = narrowest.min(smallest_margin);
+    let wanted = if feature.is_finite() && feature > 0.0 {
+        ((hi - lo) / (feature / 12.0)).ceil()
+    } else {
+        MIN_GRID_POINTS as f64
+    };
+
+    let points = if wanted.is_finite() {
+        (wanted as usize).clamp(MIN_GRID_POINTS, MAX_GRID_POINTS)
+    } else {
+        MIN_GRID_POINTS
+    };
+
+    (lo, hi, points)
+}
+
+/// Densities of each team sampled on the shared grid.
+struct Sampled {
+    lo: f64,
+    step: f64,
+    points: usize,
+    density: Vec<Vec<f64>>,
+}
+
+impl Sampled {
+    fn new(perf: &[Gaussian], margins: &Margins) -> Self {
+        let (lo, hi, points) = grid_shape(perf, margins);
+        let step = (hi - lo) / (points - 1) as f64;
+        let density = perf
+            .iter()
+            .map(|&g| {
+                (0..points)
+                    .map(|i| density(g, lo + i as f64 * step))
+                    .collect()
+            })
+            .collect();
+        Self {
+            lo,
+            step,
+            points,
+            density,
+        }
+    }
+
+    fn node(&self, i: usize) -> f64 {
+        self.lo + i as f64 * self.step
+    }
+}
+
+/// `P(order[0] >= order[1] >= ... )` with the given adjacency pattern.
+///
+/// `tied[k]` says whether `order[k]` and `order[k + 1]` finish within that
+/// pair's draw margin. The recursion runs bottom-up: `carry` holds, for each
+/// grid node, the probability that everything *below* the current team holds
+/// given that team landed on that node. A strict gap reads a cumulative
+/// integral; a tie reads a window. Both are O(1) against one prefix array,
+/// so each level costs O(grid) and the whole order costs O(teams * grid).
+fn order_probability(margins: &Margins, sampled: &Sampled, order: &[usize], tied: &[bool]) -> f64 {
+    let mut carry = vec![1.0; sampled.points];
+
+    for k in (0..order.len() - 1).rev() {
+        let below = order[k + 1];
+        let above = order[k];
+        let margin = margins.get(above, below);
+
+        let integrand: Vec<f64> = (0..sampled.points)
+            .map(|i| sampled.density[below][i] * carry[i])
+            .collect();
+        let cumulative = quadrature::Grid::from_values(sampled.lo, sampled.step, integrand);
+
+        carry = (0..sampled.points)
+            .map(|i| {
+                let x = sampled.node(i);
+                if tied[k] {
+                    // Sorted order already implies `below <= above`, so the
+                    // tie window is one-sided: [x - margin, x].
+                    cumulative.integral_between(x - margin, x)
+                } else {
+                    cumulative.integral_to(x - margin)
+                }
+            })
+            .collect();
+    }
+
+    let top = order[0];
+    let integrand: Vec<f64> = (0..sampled.points)
+        .map(|i| sampled.density[top][i] * carry[i])
+        .collect();
+    quadrature::Grid::from_values(sampled.lo, sampled.step, integrand).total()
+}
+
+/// Dense ranks implied by a sorted order and its tie pattern.
+fn ranks_of(order: &[usize], tied: &[bool], n: usize) -> Vec<u32> {
+    let mut ranks = vec![0u32; n];
+    let mut rank = 0u32;
+    ranks[order[0]] = 0;
+    for k in 0..order.len() - 1 {
+        if !tied[k] {
+            rank += 1;
+        }
+        ranks[order[k + 1]] = rank;
+    }
+    ranks
+}
+
+/// Every (order, tie-pattern) event, or only the strict ones when no pair can
+/// draw — a tie then has probability exactly zero and is not worth integrating.
+fn events(n: usize, strict_only: bool) -> Vec<(Vec<usize>, Vec<bool>)> {
+    fn permute(current: &mut Vec<usize>, k: usize, out: &mut Vec<Vec<usize>>) {
+        if k == current.len() {
+            out.push(current.clone());
+            return;
+        }
+        for i in k..current.len() {
+            current.swap(k, i);
+            permute(current, k + 1, out);
+            current.swap(k, i);
+        }
+    }
+
+    let mut orders = Vec::new();
+    permute(&mut (0..n).collect(), 0, &mut orders);
+
+    let patterns: Vec<Vec<bool>> = if strict_only {
+        vec![vec![false; n - 1]]
+    } else {
+        (0..(1u32 << (n - 1)))
+            .map(|mask| (0..n - 1).map(|i| mask >> i & 1 == 1).collect())
+            .collect()
+    };
+
+    let mut out = Vec::with_capacity(orders.len() * patterns.len());
+    for order in orders {
+        for pattern in &patterns {
+            out.push((order.clone(), pattern.clone()));
+        }
+    }
+    out
+}
+
+/// The full distribution over finishing orders, aggregated by rank vector.
+///
+/// Orders that differ only *within* a tied group describe the same finishing
+/// order, so their probabilities are summed into one entry.
+pub(crate) fn outcome_distribution(perf: &[Gaussian], margins: &Margins) -> Vec<(Vec<u32>, f64)> {
+    let n = perf.len();
+    let sampled = Sampled::new(perf, margins);
+
+    let mut aggregated: Vec<(Vec<u32>, f64)> = Vec::new();
+    for (order, tied) in events(n, margins.all_zero()) {
+        let p = order_probability(margins, &sampled, &order, &tied);
+        let ranks = ranks_of(&order, &tied, n);
+        match aggregated.iter_mut().find(|(r, _)| *r == ranks) {
+            Some((_, acc)) => *acc += p,
+            None => aggregated.push((ranks, p)),
+        }
+    }
+
+    aggregated.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap_or(std::cmp::Ordering::Equal));
+    aggregated
+}
+
+/// All permutations of `items`.
+fn permutations(items: &[usize]) -> Vec<Vec<usize>> {
+    fn go(current: &mut Vec<usize>, k: usize, out: &mut Vec<Vec<usize>>) {
+        if k == current.len() {
+            out.push(current.clone());
+            return;
+        }
+        for i in k..current.len() {
+            current.swap(k, i);
+            go(current, k + 1, out);
+            current.swap(k, i);
+        }
+    }
+    let mut out = Vec::new();
+    go(&mut items.to_vec(), 0, &mut out);
+    out
+}
+
+/// Every (order, tie-pattern) event consistent with a grouping by rank.
+///
+/// Teams sharing a rank may finish in any internal order, so this is the
+/// product of each group's permutations. Adjacencies inside a group are ties;
+/// the adjacency joining one group to the next is not.
+fn orders_for_groups(groups: &[Vec<usize>]) -> Vec<(Vec<usize>, Vec<bool>)> {
+    let per_group: Vec<Vec<Vec<usize>>> = groups.iter().map(|g| permutations(g)).collect();
+
+    let mut out = Vec::new();
+    let mut choice = vec![0usize; groups.len()];
+
+    loop {
+        let mut order = Vec::new();
+        let mut tied = Vec::new();
+        for (gi, group) in per_group.iter().enumerate() {
+            for (offset, &member) in group[choice[gi]].iter().enumerate() {
+                if !order.is_empty() {
+                    tied.push(offset != 0);
+                }
+                order.push(member);
+            }
+        }
+        out.push((order, tied));
+
+        let mut k = 0;
+        loop {
+            if k == choice.len() {
+                return out;
+            }
+            choice[k] += 1;
+            if choice[k] < per_group[k].len() {
+                break;
+            }
+            choice[k] = 0;
+            k += 1;
+        }
+    }
+}
+
+/// Probability of one specific rank vector.
+///
+/// Ties in `ranks` mean the tied teams may finish in any internal order, so
+/// this sums the orders consistent with the requested ranking rather than
+/// picking one.
+pub(crate) fn ranking_probability(perf: &[Gaussian], margins: &Margins, ranks: &[u32]) -> f64 {
+    let n = perf.len();
+    let sampled = Sampled::new(perf, margins);
+
+    let mut distinct: Vec<u32> = ranks.to_vec();
+    distinct.sort_unstable();
+    distinct.dedup();
+
+    let groups: Vec<Vec<usize>> = distinct
+        .iter()
+        .map(|&r| (0..n).filter(|&i| ranks[i] == r).collect())
+        .collect();
+
+    orders_for_groups(&groups)
+        .iter()
+        .map(|(order, tied)| order_probability(margins, &sampled, order, tied))
+        .sum()
+}
+
+/// A distribution over the ways a contest could finish.
+///
+/// Each entry pairs a rank vector — the same shape [`crate::Outcome::ranking`]
+/// takes, with equal ranks meaning a tie — against its probability. Entries
+/// are ordered most likely first, and cover the whole outcome space, so the
+/// probabilities sum to one.
+///
+/// The rank vectors compose directly with inference: feeding one to
+/// `Game::ranked` asks "what would we believe if *this* happened", which is
+/// what an expected-information-gain calculation needs alongside the weight.
+#[derive(Clone, Debug, PartialEq)]
+pub struct Prediction {
+    outcomes: Vec<(Vec<u32>, f64)>,
+}
+
+impl Prediction {
+    pub(crate) fn new(outcomes: Vec<(Vec<u32>, f64)>) -> Self {
+        Self { outcomes }
+    }
+
+    /// Every possible finishing order and its probability, most likely first.
+    pub fn outcomes(&self) -> impl ExactSizeIterator<Item = (&[u32], f64)> {
+        self.outcomes.iter().map(|(r, p)| (r.as_slice(), *p))
+    }
+
+    /// The single most likely finishing order.
+    #[must_use]
+    pub fn most_likely(&self) -> Option<(&[u32], f64)> {
+        self.outcomes.first().map(|(r, p)| (r.as_slice(), *p))
+    }
+
+    /// Probability of one specific finishing order, or zero if it cannot occur.
+    #[must_use]
+    pub fn probability_of(&self, ranks: &[u32]) -> f64 {
+        self.outcomes
+            .iter()
+            .find(|(r, _)| r.as_slice() == ranks)
+            .map_or(0.0, |(_, p)| *p)
+    }
+
+    /// `P(team i finishes strictly first)`, for each team.
+    ///
+    /// Sums to less than one exactly when the top place can be shared; the
+    /// shortfall is [`Prediction::shared_first_place`].
+    #[must_use]
+    pub fn win_probabilities(&self) -> Vec<f64> {
+        let n = self.outcomes.first().map_or(0, |(r, _)| r.len());
+        let mut wins = vec![0.0; n];
+        for (ranks, p) in &self.outcomes {
+            let leaders = ranks.iter().filter(|&&r| r == 0).count();
+            if leaders == 1 {
+                let winner = ranks.iter().position(|&r| r == 0).expect("a rank-0 team");
+                wins[winner] += p;
+            }
+        }
+        wins
+    }
+
+    /// Probability that two or more teams share first place.
+    #[must_use]
+    pub fn shared_first_place(&self) -> f64 {
+        self.outcomes
+            .iter()
+            .filter(|(r, _)| r.iter().filter(|&&x| x == 0).count() > 1)
+            .map(|(_, p)| p)
+            .sum()
+    }
+
+    /// Total probability mass, which should be one.
+    ///
+    /// Exposed because it is a genuine check on the numerics rather than a
+    /// formality: the outcome space is exhaustive and disjoint by construction,
+    /// so any drift from one is integration error and nothing else.
+    #[must_use]
+    pub fn total(&self) -> f64 {
+        self.outcomes.iter().map(|(_, p)| p).sum()
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    fn g(mu: f64, sigma: f64) -> Gaussian {
+        Gaussian::from_ms(mu, sigma)
+    }
+
+    fn flat(n: usize, eps: f64) -> Margins {
+        Margins::new(n, |_, _| eps)
+    }
+
+    /// Exact two-team result: `P(a first) = Phi((mu_a - mu_b - eps) / sd)`.
+    fn closed_form_two(a: Gaussian, b: Gaussian, eps: f64) -> (f64, f64) {
+        let sd = (a.sigma().powi(2) + b.sigma().powi(2)).sqrt();
+        (
+            phi((a.mu() - b.mu() - eps) / sd),
+            phi((b.mu() - a.mu() - eps) / sd),
+        )
+    }
+
+    #[test]
+    fn two_team_win_probabilities_match_the_closed_form() {
+        for (ma, sa, mb, sb, eps) in [
+            (0.0, 6.0, 0.0, 6.0, 0.0),
+            (3.0, 6.0, -2.0, 1.0, 0.0),
+            (0.0, 6.0, 0.0, 6.0, 2.0),
+            (3.0, 6.0, -2.0, 1.0, 1.5),
+            (40.0, 1.0, 0.0, 1.0, 0.0),
+        ] {
+            let perf = [g(ma, sa), g(mb, sb)];
+            let got = win_probabilities(&perf, &flat(2, eps));
+            let (wa, wb) = closed_form_two(perf[0], perf[1], eps);
+            assert!(
+                (got[0] - wa).abs() < 1e-7 && (got[1] - wb).abs() < 1e-7,
+                "mu=({ma},{mb}) sigma=({sa},{sb}) eps={eps}: got {got:?}, want [{wa}, {wb}]"
+            );
+        }
+    }
+
+    /// The identity that a wrong-but-plausible implementation cannot fake:
+    /// with no draw margin, exactly one team finishes first.
+    #[test]
+    fn win_probabilities_sum_to_one_without_a_draw_margin() {
+        for perf in [
+            vec![g(0.0, 6.0), g(0.0, 6.0)],
+            vec![g(5.0, 6.0), g(0.0, 3.0), g(-5.0, 1.0)],
+            vec![
+                g(8.0, 2.0),
+                g(3.0, 6.0),
+                g(0.0, 1.0),
+                g(-3.0, 4.0),
+                g(-8.0, 6.0),
+            ],
+        ] {
+            let sum: f64 = win_probabilities(&perf, &flat(perf.len(), 0.0))
+                .iter()
+                .sum();
+            assert!(
+                (sum - 1.0).abs() < 1e-7,
+                "{} teams: sum = {sum}",
+                perf.len()
+            );
+        }
+    }
+
+    /// A rival with a tiny sigma is a step function in disguise. Fixed-node
+    /// quadrature steps over it and lands ~1e-2 out while still looking like a
+    /// probability; this is the case that rules that approach out.
+    #[test]
+    fn win_probabilities_survive_a_rival_with_a_tiny_sigma() {
+        let perf = [g(0.0, 0.001), g(0.5, 6.0), g(-0.5, 6.0)];
+        let got = win_probabilities(&perf, &flat(3, 0.0));
+        let sum: f64 = got.iter().sum();
+        assert!((sum - 1.0).abs() < 1e-6, "sum = {sum}, probs = {got:?}");
+    }
+
+    #[test]
+    fn a_stronger_team_is_more_likely_to_win() {
+        let perf = [g(10.0, 3.0), g(0.0, 3.0), g(-10.0, 3.0)];
+        let p = win_probabilities(&perf, &flat(3, 0.0));
+        assert!(p[0] > p[1] && p[1] > p[2], "not monotone: {p:?}");
+    }
+
+    #[test]
+    fn identical_teams_are_equally_likely_to_win() {
+        let perf = [g(1.0, 4.0), g(1.0, 4.0), g(1.0, 4.0)];
+        let p = win_probabilities(&perf, &flat(3, 0.0));
+        for probs in p.windows(2) {
+            assert!((probs[0] - probs[1]).abs() < 1e-9, "asymmetric: {p:?}");
+        }
+    }
+
+    /// Every realisation sorts into exactly one finishing order, so the whole
+    /// distribution must sum to one — with or without a draw margin.
+    #[test]
+    fn outcome_distribution_sums_to_one() {
+        for (perf, eps) in [
+            (vec![g(0.0, 6.0), g(0.0, 6.0)], 0.0),
+            (vec![g(0.0, 6.0), g(0.0, 6.0)], 2.0),
+            (vec![g(0.0, 6.0), g(0.0, 6.0), g(0.0, 6.0)], 0.0),
+            (vec![g(5.0, 6.0), g(0.0, 3.0), g(-5.0, 1.0)], 1.5),
+            (vec![g(0.0, 0.05), g(0.5, 6.0), g(-0.5, 6.0)], 1.0),
+            (
+                vec![g(6.0, 2.0), g(2.0, 6.0), g(-2.0, 1.0), g(-6.0, 4.0)],
+                1.0,
+            ),
+        ] {
+            let n = perf.len();
+            let dist = outcome_distribution(&perf, &flat(n, eps));
+            let sum: f64 = dist.iter().map(|(_, p)| p).sum();
+            assert!(
+                (sum - 1.0).abs() < 1e-6,
+                "{n} teams, eps={eps}: sum = {sum} over {} outcomes",
+                dist.len()
+            );
+            assert!(dist.iter().all(|(_, p)| *p >= 0.0), "negative probability");
+        }
+    }
+
+    /// With two teams the distribution is the exact win/draw/loss triple.
+    #[test]
+    fn two_team_distribution_matches_the_closed_form() {
+        let perf = [g(3.0, 6.0), g(-2.0, 1.0)];
+        let eps = 1.5;
+        let dist = outcome_distribution(&perf, &flat(2, eps));
+        let (wa, wb) = closed_form_two(perf[0], perf[1], eps);
+
+        let find = |ranks: &[u32]| {
+            dist.iter()
+                .find(|(r, _)| r == ranks)
+                .map_or(0.0, |(_, p)| *p)
+        };
+
+        assert!(
+            (find(&[0, 1]) - wa).abs() < 1e-6,
+            "a wins: {}",
+            find(&[0, 1])
+        );
+        assert!(
+            (find(&[1, 0]) - wb).abs() < 1e-6,
+            "b wins: {}",
+            find(&[1, 0])
+        );
+        assert!(
+            (find(&[0, 0]) - (1.0 - wa - wb)).abs() < 1e-6,
+            "draw: {}",
+            find(&[0, 0])
+        );
+    }
+
+    /// Asking for one ranking must agree with that ranking's entry in the
+    /// full distribution — the two use different code paths to the same value.
+    #[test]
+    fn ranking_probability_agrees_with_the_distribution() {
+        let perf = [g(5.0, 6.0), g(0.0, 3.0), g(-5.0, 1.0)];
+        let eps = 1.5;
+        let margins = flat(3, eps);
+        let dist = outcome_distribution(&perf, &margins);
+
+        for (ranks, expected) in &dist {
+            let direct = ranking_probability(&perf, &margins, ranks);
+            assert!(
+                (direct - expected).abs() < 1e-9,
+                "ranks {ranks:?}: direct {direct} vs distribution {expected}"
+            );
+        }
+    }
+
+    /// Tie mass is controlled by the draw margin. Only the *all-tied* outcome
+    /// is monotone in it: every one of its constraints is a window that widens
+    /// with the margin. A partially-tied outcome like `[0, 0, 1]` is not, and
+    /// must not be asserted to be — widening the margin makes its tie easier
+    /// but its "and the last team is strictly behind by more than the margin"
+    /// clause harder, so it peaks and then falls.
+    #[test]
+    fn all_tied_probability_grows_with_the_draw_margin() {
+        let perf = [g(0.0, 4.0), g(0.0, 4.0), g(-8.0, 2.0)];
+        let mut previous = 0.0;
+        for eps in [0.0, 0.5, 1.0, 2.0, 4.0, 8.0, 24.0] {
+            let p = ranking_probability(&perf, &flat(3, eps), &[0, 0, 0]);
+            assert!(p >= previous, "eps={eps}: {p} < {previous}");
+            if eps == 0.0 {
+                assert!(p < 1e-12, "a tie needs a margin, got {p}");
+            }
+            previous = p;
+        }
+        assert!(
+            previous > 0.9,
+            "a very wide margin ties everyone: {previous}"
+        );
+    }
+
+    /// The converse, stated as the non-property it is: a partially-tied
+    /// outcome is non-monotone in the margin. Pinning this down stops a future
+    /// change from "fixing" it into monotonicity and quietly breaking the model.
+    #[test]
+    fn a_partially_tied_outcome_peaks_in_the_middle() {
+        let perf = [g(0.0, 4.0), g(0.0, 4.0), g(-8.0, 2.0)];
+        let sweep: Vec<f64> = [0.5, 2.0, 4.0, 8.0, 16.0]
+            .iter()
+            .map(|&eps| ranking_probability(&perf, &flat(3, eps), &[0, 0, 1]))
+            .collect();
+        let peak = sweep
+            .iter()
+            .enumerate()
+            .fold(
+                (0, 0.0),
+                |(bi, bv), (i, &v)| if v > bv { (i, v) } else { (bi, bv) },
+            )
+            .0;
+        assert!(
+            peak > 0 && peak < sweep.len() - 1,
+            "expected an interior peak: {sweep:?}"
+        );
+    }
+
+    /// With no draw margin a tie has probability exactly zero, and the
+    /// enumeration must not waste work pretending otherwise.
+    #[test]
+    fn ties_are_impossible_without_a_draw_margin() {
+        let perf = [g(0.0, 4.0), g(0.0, 4.0), g(0.0, 4.0)];
+        let dist = outcome_distribution(&perf, &flat(3, 0.0));
+        assert_eq!(dist.len(), 6, "expected only the 6 strict orders: {dist:?}");
+        assert!(dist.iter().all(|(r, _)| {
+            let mut seen = r.clone();
+            seen.sort_unstable();
+            seen.dedup();
+            seen.len() == r.len()
+        }));
+    }
+}

src/quadrature.rs

@@ -0,0 +1,322 @@
+//! Deterministic numerical integration for the prediction paths.
+//!
+//! Prediction asks two questions that have no closed form beyond two teams:
+//! "who finishes first" and "how likely is this exact finishing order". Both
+//! reduce to integrals over a single performance variable, so neither needs a
+//! sampler — and that matters, because a Monte Carlo predictor would make
+//! `predict_*` non-reproducible and would answer a slightly different question
+//! on every call.
+//!
+//! Two routines live here:
+//!
+//! - [`integrate`], adaptive Gauss-Kronrod G7-K15, for the first-place
+//!   marginals. It carries its own error estimate, so it can refine where the
+//!   integrand actually bends instead of guessing a node count up front.
+//! - [`Grid`], a uniform grid with trapezoid prefix sums, for the ranking
+//!   chain recursion, where each level needs the *running* integral of the
+//!   level below at arbitrary points rather than one definite integral.
+//!
+//! Fixed-node Gauss-Hermite is the obvious tool for the first of these and is
+//! a trap: the integrand is a product of normal CDFs, and when one team's
+//! sigma is much smaller than the integrating team's, that product turns into
+//! a near-step function narrower than the node spacing. The nodes step over
+//! it and the result is wrong by ~1e-2 while still looking like a probability.
+//! Adaptive refinement is what makes the small-sigma case safe.
+
+/// Kronrod 15-point abscissae, non-negative half, descending.
+const XGK: [f64; 8] = [
+    0.991_455_371_120_813,
+    0.949_107_912_342_759,
+    0.864_864_423_359_769,
+    0.741_531_185_599_394,
+    0.586_087_235_467_691,
+    0.405_845_151_377_397,
+    0.207_784_955_007_898,
+    0.0,
+];
+
+/// Kronrod 15-point weights, matching [`XGK`].
+const WGK: [f64; 8] = [
+    0.022_935_322_010_529,
+    0.063_092_092_629_979,
+    0.104_790_010_322_250,
+    0.140_653_259_715_525,
+    0.169_004_726_639_267,
+    0.190_350_578_064_785,
+    0.204_432_940_075_298,
+    0.209_482_141_084_728,
+];
+
+/// Gauss 7-point weights, applying to the odd-indexed [`XGK`] entries.
+const WG: [f64; 4] = [
+    0.129_484_966_168_870,
+    0.279_705_391_489_277,
+    0.381_830_050_505_119,
+    0.417_959_183_673_469,
+];
+
+/// Panels are bisected worst-first; this bounds the work on a pathological
+/// integrand rather than letting it spin.
+const MAX_SUBDIVISIONS: usize = 200;
+
+/// One G7-K15 panel over `[a, b]`: `(integral, absolute error estimate)`.
+///
+/// The error estimate is the gap between the embedded 7-point Gauss rule and
+/// the 15-point Kronrod extension. It is the only reason this is preferable
+/// to a fixed rule: it tells the caller *where* the integrand is hard.
+fn gk15<F: Fn(f64) -> f64>(f: &F, a: f64, b: f64) -> (f64, f64) {
+    let centre = 0.5 * (a + b);
+    let half = 0.5 * (b - a);
+
+    let mut kronrod = 0.0;
+    let mut gauss = 0.0;
+
+    for i in 0..8 {
+        let offset = XGK[i] * half;
+        // XGK[7] is the centre node and must not be counted twice.
+        let sum = if i == 7 {
+            f(centre)
+        } else {
+            f(centre - offset) + f(centre + offset)
+        };
+        kronrod += WGK[i] * sum;
+        if i % 2 == 1 {
+            gauss += WG[i / 2] * sum;
+        }
+    }
+
+    (kronrod * half, ((kronrod - gauss) * half).abs())
+}
+
+/// Adaptively integrate `f` over `[a, b]` to relative tolerance `tol`.
+///
+/// `seeds` are interior points where the integrand is known to bend sharply —
+/// for a product of normal CDFs, each rival's transition centre. Splitting
+/// there up front costs nothing and saves the adaptive loop from having to
+/// discover a step by bisection.
+///
+/// Returns the integral. The error estimate is consumed internally rather
+/// than returned: callers here integrate probability densities, where the
+/// meaningful check is the sum-to-one identity over a whole outcome space,
+/// not a per-integral residual.
+pub(crate) fn integrate<F: Fn(f64) -> f64>(f: F, a: f64, b: f64, seeds: &[f64], tol: f64) -> f64 {
+    // Explicit rather than `!(b > a)`: a NaN bound must fall through to zero
+    // rather than being read as a valid ordering.
+    if a.partial_cmp(&b) != Some(std::cmp::Ordering::Less) {
+        return 0.0;
+    }
+
+    let mut edges: Vec<f64> = Vec::with_capacity(seeds.len() + 2);
+    edges.push(a);
+    edges.push(b);
+    for &s in seeds {
+        if s > a && s < b {
+            edges.push(s);
+        }
+    }
+    edges.sort_by(|p, q| p.partial_cmp(q).expect("integration bounds are finite"));
+    edges.dedup();
+
+    // (lo, hi, integral, error)
+    let mut panels: Vec<(f64, f64, f64, f64)> = edges
+        .windows(2)
+        .map(|w| {
+            let (v, e) = gk15(&f, w[0], w[1]);
+            (w[0], w[1], v, e)
+        })
+        .collect();
+
+    for _ in 0..MAX_SUBDIVISIONS {
+        let total: f64 = panels.iter().map(|p| p.2).sum();
+        let error: f64 = panels.iter().map(|p| p.3).sum();
+
+        // Absolute floor as well as relative: these integrands are
+        // probabilities, so an absolute 1e-15 is already past the useful
+        // precision of the underlying `cdf`.
+        if error <= tol * total.abs().max(1e-12) || error < 1e-15 {
+            break;
+        }
+
+        let worst = panels
+            .iter()
+            .enumerate()
+            .fold((0usize, f64::NEG_INFINITY), |(bi, be), (i, p)| {
+                if p.3 > be { (i, p.3) } else { (bi, be) }
+            })
+            .0;
+
+        let (lo, hi, _, _) = panels[worst];
+        let mid = 0.5 * (lo + hi);
+        // Bisection has hit the floating-point floor; refining further would
+        // loop without reducing the error.
+        if !(mid > lo && mid < hi) {
+            break;
+        }
+
+        let (v1, e1) = gk15(&f, lo, mid);
+        let (v2, e2) = gk15(&f, mid, hi);
+        panels[worst] = (lo, mid, v1, e1);
+        panels.push((mid, hi, v2, e2));
+    }
+
+    panels.iter().map(|p| p.2).sum()
+}
+
+/// A uniform grid carrying trapezoid prefix sums of one integrand.
+///
+/// The ranking recursion needs, at every level, the running integral of the
+/// level below evaluated at arbitrary points — a cumulative integral, not a
+/// definite one. Prefix sums give that in O(1) per query after an O(G) build,
+/// which is what keeps a full ranking probability linear in the team count.
+pub(crate) struct Grid {
+    lo: f64,
+    step: f64,
+    /// Integrand sampled at each node.
+    values: Vec<f64>,
+    /// `prefix[i]` is the integral from `lo` to node `i`.
+    prefix: Vec<f64>,
+}
+
+impl Grid {
+    /// Build directly from already-sampled values.
+    ///
+    /// The ranking recursion evaluates every level on the same nodes, so the
+    /// per-team densities are sampled once and reused; re-evaluating `exp`
+    /// per level would dominate the cost.
+    pub(crate) fn from_values(lo: f64, step: f64, values: Vec<f64>) -> Self {
+        let mut prefix = vec![0.0; values.len()];
+        for i in 1..values.len() {
+            prefix[i] = prefix[i - 1] + 0.5 * step * (values[i - 1] + values[i]);
+        }
+        Self {
+            lo,
+            step,
+            values,
+            prefix,
+        }
+    }
+
+    /// Integral from the grid's lower bound up to `x`.
+    ///
+    /// Clamped at both ends: the caller sizes the grid to cover the whole
+    /// support, so a query outside it is asking for a tail that is zero (below)
+    /// or the whole mass (above).
+    pub(crate) fn integral_to(&self, x: f64) -> f64 {
+        let last = self.values.len() - 1;
+        if x <= self.lo {
+            return 0.0;
+        }
+        if x >= self.lo + last as f64 * self.step {
+            return self.prefix[last];
+        }
+
+        let scaled = (x - self.lo) / self.step;
+        let i = scaled.floor() as usize;
+        let frac = scaled - i as f64;
+
+        // Whole cells, plus the trapezoid over the partial cell. The integrand
+        // is linear within a cell under the trapezoid rule, so the partial
+        // piece is exact with respect to that same approximation.
+        self.prefix[i]
+            + frac
+                * self.step
+                * (self.values[i] + 0.5 * frac * (self.values[i + 1] - self.values[i]))
+    }
+
+    /// Integral over `[from, to]`.
+    pub(crate) fn integral_between(&self, from: f64, to: f64) -> f64 {
+        (self.integral_to(to) - self.integral_to(from)).max(0.0)
+    }
+
+    /// Total integral over the whole grid.
+    pub(crate) fn total(&self) -> f64 {
+        self.prefix[self.values.len() - 1]
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    const TOL: f64 = 1e-10;
+
+    /// Sample `f` over `[lo, hi]` at `points` nodes.
+    fn sample<F: FnMut(f64) -> f64>(lo: f64, hi: f64, points: usize, mut f: F) -> Grid {
+        let step = (hi - lo) / (points - 1) as f64;
+        Grid::from_values(
+            lo,
+            step,
+            (0..points).map(|i| f(lo + i as f64 * step)).collect(),
+        )
+    }
+
+    #[test]
+    fn integrates_a_polynomial_exactly() {
+        // G7-K15 is exact for polynomials well past cubic, so a single panel
+        // should already be at round-off.
+        let v = integrate(|x| 3.0 * x * x + 2.0 * x + 1.0, 0.0, 2.0, &[], TOL);
+        assert!((v - 14.0).abs() < 1e-12, "got {v}");
+    }
+
+    #[test]
+    fn integrates_a_gaussian_density_to_one() {
+        let f = |x: f64| (-0.5 * x * x).exp() / (2.0 * std::f64::consts::PI).sqrt();
+        let v = integrate(f, -10.0, 10.0, &[], TOL);
+        assert!((v - 1.0).abs() < 1e-12, "got {v}");
+    }
+
+    #[test]
+    fn resolves_a_step_far_narrower_than_the_initial_panel() {
+        // The failure mode that rules out fixed-node quadrature: a transition
+        // 1e-4 wide inside a range of 20. A fixed rule steps over it.
+        let f = |x: f64| if x < 0.5 { 0.0 } else { 1.0 };
+        let v = integrate(f, -10.0, 10.0, &[0.5], TOL);
+        assert!((v - 9.5).abs() < 1e-6, "got {v}");
+    }
+
+    #[test]
+    fn seeds_do_not_change_the_value_of_a_smooth_integrand() {
+        let f = |x: f64| (-0.5 * x * x).exp();
+        let plain = integrate(f, -8.0, 8.0, &[], TOL);
+        let seeded = integrate(f, -8.0, 8.0, &[-3.0, 0.25, 5.5], TOL);
+        assert!((plain - seeded).abs() < 1e-12, "{plain} vs {seeded}");
+    }
+
+    #[test]
+    fn empty_or_inverted_range_integrates_to_zero() {
+        assert_eq!(integrate(|_| 1.0, 1.0, 1.0, &[], TOL), 0.0);
+        assert_eq!(integrate(|_| 1.0, 2.0, 1.0, &[], TOL), 0.0);
+    }
+
+    #[test]
+    fn grid_prefix_matches_a_known_cumulative_integral() {
+        // f(x) = x over [0, 4]; integral to x is x^2/2.
+        let g = sample(0.0, 4.0, 4001, |x| x);
+        for probe in [0.0, 0.5, 1.0, 2.5, 3.75, 4.0] {
+            let want = probe * probe / 2.0;
+            let got = g.integral_to(probe);
+            assert!(
+                (got - want).abs() < 1e-9,
+                "at {probe}: got {got}, want {want}"
+            );
+        }
+        assert!((g.total() - 8.0).abs() < 1e-9);
+    }
+
+    #[test]
+    fn grid_between_is_the_difference_of_two_prefixes() {
+        let g = sample(-5.0, 5.0, 8001, |x| (-0.5 * x * x).exp());
+        let whole = g.integral_between(-5.0, 5.0);
+        let split = g.integral_between(-5.0, 0.3) + g.integral_between(0.3, 5.0);
+        assert!((whole - split).abs() < 1e-12, "{whole} vs {split}");
+    }
+
+    #[test]
+    fn grid_clamps_queries_outside_its_support() {
+        let g = sample(0.0, 1.0, 101, |_| 1.0);
+        assert_eq!(g.integral_to(-3.0), 0.0);
+        assert!((g.integral_to(9.0) - 1.0).abs() < 1e-12);
+        // Reversed bounds must not produce negative probability mass.
+        assert_eq!(g.integral_between(0.8, 0.2), 0.0);
+    }
+}

tests/api_shape.rs

@@ -203,7 +203,7 @@ fn predict_quality_two_teams() {
     h.record_winner(&"a", &"b", 1).unwrap();
     h.converge().unwrap();
 
-    let q = h.predict_quality(&[&[&"a"], &[&"b"]]);
+    let q = h.predict_quality(&[&[&"a"], &[&"b"]]).unwrap();
     assert!(q > 0.0 && q <= 1.0);
 }
 
@@ -219,10 +219,14 @@ fn predict_outcome_two_teams_sums_to_one() {
     h.record_winner(&"a", &"b", 1).unwrap();
     h.converge().unwrap();
 
-    let p = h.predict_outcome(&[&[&"a"], &[&"b"]]);
+    let p = h.predict_outcome(&[&[&"a"], &[&"b"]]).unwrap();
-    assert_eq!(p.len(), 2);
+    let wins = p.win_probabilities();
-    assert!((p[0] + p[1] - 1.0).abs() < 1e-9);
+    assert_eq!(wins.len(), 2);
-    assert!(p[0] > p[1]);
+    // With p_draw == 0 there is no draw outcome, so the two win
+    // probabilities are the whole space.
+    assert!((p.total() - 1.0).abs() < 1e-9, "total = {}", p.total());
+    assert!((wins[0] + wins[1] - 1.0).abs() < 1e-9);
+    assert!(wins[0] > wins[1]);
 }
 
 #[test]

tests/prediction.rs

@@ -0,0 +1,214 @@
+//! Prediction API: N-team outcomes, draw mass, and the error paths that used
+//! to be panics or silent wrong answers.
+
+use trueskill_tt::{History, InferenceError, MAX_PREDICTED_TEAMS};
+
+fn history_with(names: &[&'static str], p_draw: f64) -> History {
+    let mut h = History::builder().p_draw(p_draw).build();
+    // Give every competitor a recorded skill by playing a small round robin.
+    for pair in names.windows(2) {
+        h.record_winner(&pair[0], &pair[1], 1).unwrap();
+    }
+    h.converge().unwrap();
+    h
+}
+
+#[test]
+fn unknown_keys_are_reported_not_silently_dropped() {
+    let h = history_with(&["a", "b"], 0.0);
+
+    let err = h
+        .predict_outcome(&[&[&"a"], &[&"ghost"]])
+        .expect_err("an unknown key must not yield a confident prediction");
+    assert_eq!(err, InferenceError::UnknownKey { team: 1, member: 0 });
+
+    // Every prediction entry point, not just one.
+    assert!(
+        h.predict_win_probabilities(&[&[&"a"], &[&"ghost"]])
+            .is_err()
+    );
+    assert!(h.predict_quality(&[&[&"a"], &[&"ghost"]]).is_err());
+    assert!(h.predict_ranking(&[&[&"a"], &[&"ghost"]], &[0, 1]).is_err());
+}
+
+#[test]
+fn an_entirely_unknown_team_is_an_error() {
+    let h = history_with(&["a", "b"], 0.0);
+    let err = h.predict_outcome(&[&[&"a"], &[&"x", &"y"]]).unwrap_err();
+    assert_eq!(err, InferenceError::UnknownKey { team: 1, member: 0 });
+}
+
+#[test]
+fn degenerate_team_shapes_are_errors_rather_than_panics() {
+    let h = history_with(&["a", "b"], 0.0);
+
+    assert_eq!(
+        h.predict_outcome(&[&[&"a"]]).unwrap_err(),
+        InferenceError::NotEnoughTeams { got: 1 }
+    );
+    assert_eq!(
+        h.predict_outcome(&[]).unwrap_err(),
+        InferenceError::NotEnoughTeams { got: 0 }
+    );
+    assert_eq!(
+        h.predict_outcome(&[&[&"a"], &[]]).unwrap_err(),
+        InferenceError::EmptyTeam { team: 1 }
+    );
+}
+
+#[test]
+fn more_than_two_teams_no_longer_panics() {
+    let h = history_with(&["a", "b", "c"], 0.0);
+    let p = h
+        .predict_outcome(&[&[&"a"], &[&"b"], &[&"c"]])
+        .expect("three teams must be supported");
+    assert!((p.total() - 1.0).abs() < 1e-6, "total = {}", p.total());
+    // Three teams, no draws possible: exactly the six strict orderings.
+    assert_eq!(p.outcomes().len(), 6);
+}
+
+#[test]
+fn the_outcome_space_is_capped_rather_than_hanging() {
+    let names: Vec<&'static str> = vec!["a", "b", "c", "d", "e", "f", "g", "h"];
+    let h = history_with(&names, 0.0);
+
+    let teams: Vec<&[&&'static str]> = Vec::new();
+    let _ = teams;
+
+    let too_many: Vec<Vec<&&str>> = names.iter().map(|n| vec![n]).collect();
+    let refs: Vec<&[&&str]> = too_many.iter().map(Vec::as_slice).collect();
+
+    let err = h.predict_outcome(&refs).unwrap_err();
+    assert_eq!(
+        err,
+        InferenceError::TooManyTeams {
+            got: 8,
+            max: MAX_PREDICTED_TEAMS
+        }
+    );
+
+    // The cheap paths stay available at any size.
+    let wins = h.predict_win_probabilities(&refs).unwrap();
+    assert_eq!(wins.len(), 8);
+    assert!(
+        (wins.iter().sum::<f64>() - 1.0).abs() < 1e-6,
+        "win probabilities must still sum to one: {wins:?}"
+    );
+}
+
+/// The defect that made every draw-enabled prediction wrong: `[p, 1 - p]`
+/// allocated no mass to a draw even with `p_draw > 0`.
+#[test]
+fn a_draw_carries_probability_mass_when_p_draw_is_positive() {
+    let h = history_with(&["a", "b"], 0.25);
+    let p = h.predict_outcome(&[&[&"a"], &[&"b"]]).unwrap();
+
+    let draw = p.probability_of(&[0, 0]);
+    assert!(draw > 0.0, "a draw-enabled model must give draws mass");
+    assert!((p.total() - 1.0).abs() < 1e-6, "total = {}", p.total());
+
+    let wins = p.win_probabilities();
+    assert!(
+        (wins.iter().sum::<f64>() + draw - 1.0).abs() < 1e-6,
+        "wins {wins:?} plus draw {draw} must be the whole space"
+    );
+    assert!(
+        (p.shared_first_place() - draw).abs() < 1e-12,
+        "a two-team draw is a shared first place"
+    );
+}
+
+#[test]
+fn a_zero_draw_probability_admits_no_ties() {
+    let h = history_with(&["a", "b"], 0.0);
+    let p = h.predict_outcome(&[&[&"a"], &[&"b"]]).unwrap();
+    assert_eq!(p.probability_of(&[0, 0]), 0.0);
+    assert!(p.shared_first_place() < 1e-12);
+}
+
+/// The two routes to a win probability run through entirely different
+/// algorithms — adaptive quadrature versus the enumerated chain recursion —
+/// so agreement between them is a real cross-check, not a tautology.
+#[test]
+fn the_cheap_and_exhaustive_paths_agree() {
+    for p_draw in [0.0, 0.1] {
+        let h = history_with(&["a", "b", "c"], p_draw);
+        let teams: &[&[&&str]] = &[&[&"a"], &[&"b"], &[&"c"]];
+
+        let cheap = h.predict_win_probabilities(teams).unwrap();
+        let exhaustive = h.predict_outcome(teams).unwrap().win_probabilities();
+
+        for (i, (a, b)) in cheap.iter().zip(&exhaustive).enumerate() {
+            assert!(
+                (a - b).abs() < 1e-6,
+                "p_draw={p_draw} team {i}: quadrature {a} vs enumeration {b}"
+            );
+        }
+    }
+}
+
+#[test]
+fn predict_ranking_agrees_with_the_distribution() {
+    let h = history_with(&["a", "b", "c"], 0.1);
+    let teams: &[&[&&str]] = &[&[&"a"], &[&"b"], &[&"c"]];
+    let dist = h.predict_outcome(teams).unwrap();
+
+    for (ranks, expected) in dist.outcomes() {
+        let direct = h.predict_ranking(teams, ranks).unwrap();
+        assert!(
+            (direct - expected).abs() < 1e-9,
+            "ranks {ranks:?}: {direct} vs {expected}"
+        );
+    }
+}
+
+#[test]
+fn predict_ranking_checks_its_shape() {
+    let h = history_with(&["a", "b"], 0.0);
+    let err = h
+        .predict_ranking(&[&[&"a"], &[&"b"]], &[0, 1, 2])
+        .unwrap_err();
+    assert!(matches!(
+        err,
+        InferenceError::MismatchedShape {
+            expected: 2,
+            got: 3,
+            ..
+        }
+    ));
+}
+
+#[test]
+fn the_stronger_competitor_is_favoured() {
+    let mut h = History::builder().build();
+    for t in 1..=10 {
+        h.record_winner(&"strong", &"weak", t).unwrap();
+    }
+    h.converge().unwrap();
+
+    let p = h.predict_outcome(&[&[&"strong"], &[&"weak"]]).unwrap();
+    let (best, _) = p.most_likely().expect("a most likely outcome");
+    assert_eq!(best, &[0, 1], "the winner should be favoured");
+
+    let wins = p.win_probabilities();
+    assert!(wins[0] > wins[1], "{wins:?}");
+}
+
+/// Unequal team sizes change the draw margin, because inference derives it
+/// from the teams' betas. Prediction has to follow, or it describes a
+/// different model than the one that will be fitted.
+#[test]
+fn team_size_affects_the_prediction() {
+    let mut h = History::builder().p_draw(0.2).build();
+    h.event(1)
+        .team(["a", "b"])
+        .team(["c"])
+        .winner(0)
+        .commit()
+        .unwrap();
+    h.converge().unwrap();
+
+    let p = h.predict_outcome(&[&[&"a", &"b"], &[&"c"]]).unwrap();
+    assert!((p.total() - 1.0).abs() < 1e-6, "total = {}", p.total());
+    assert!(p.probability_of(&[0, 0]) > 0.0);
+}

tests/quality.rs

@@ -110,7 +110,7 @@ fn history_predict_quality_supports_three_teams() {
     h.record_winner(&"b", &"c", 2).unwrap();
     h.converge().unwrap();
 
-    let q = h.predict_quality(&[&[&"a"], &[&"b"], &[&"c"]]);
+    let q = h.predict_quality(&[&[&"a"], &[&"b"], &[&"c"]]).unwrap();
     assert!(
         q.is_finite(),
         "3-team predict_quality must be finite, got {q}"