logaritmiskandClaude Opus 5
92d690d0f8
`lib.rs` advertised `Game` in "Core types" as "one match in isolation". It had no public constructor: every `Game::*` returned `OwnedGame`, so `let g: Game = Game::ranked(..)?` did not compile. Names swapped. The public type is the owned one — `Game<T, D>`, no lifetime — and the borrowing form is `pub(crate) GameRef<'a, T, D>`, which is what it always was: an implementation detail about whether the result and weight slices are borrowed from `History`'s storage. That distinction meant nothing to someone scoring one match, and it showed the module's surface twice in rustdoc, since both types carried `posteriors()` / `log_evidence()`. `one_v_one` returned `(Gaussian, Gaussian)` while every sibling returned a game, making it the one constructor you could not ask for `log_evidence()`. It returns `Self` now; `.posteriors()` recovers the old shape, and the test that covers it now also asserts the evidence of two identical ratings is exactly `ln(0.5)`. `ranked`, `scored` and `free_for_all` had `# Errors` as their entire doc, so rustdoc's index rendered the error list as the summary. They have summary lines, and `Game` has a worked example — it was advertised as a core type with none anywhere in the crate. Closes #69. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_011hcFjNDmHXZF8URGLku5zZ
1363 lines
53 KiB
Rust
1363 lines
53 KiB
Rust
//! `TrueSkill` Through Time — Bayesian skill rating over a time axis.
|
|
//!
|
|
//! Where plain `TrueSkill` gives each competitor one running estimate, `TrueSkill`
|
|
//! Through Time treats a whole history as a single model and infers skill *at
|
|
//! every point in time*. Evidence flows both directions: a result today
|
|
//! sharpens the estimate of who someone was last year, so early estimates stop
|
|
//! being frozen guesses and comparisons across eras become meaningful.
|
|
//!
|
|
//! This is a Rust port of
|
|
//! [TrueSkillThroughTime.py](https://github.com/glandfried/TrueSkillThroughTime.py).
|
|
//!
|
|
//! # Getting started
|
|
//!
|
|
//! Record results, converge, then read off skills:
|
|
//!
|
|
//! ```
|
|
//! use trueskill_tt::History;
|
|
//!
|
|
//! let mut history = History::default();
|
|
//!
|
|
//! history.record_winner(&"alice", &"bob", 1)?;
|
|
//! history.record_winner(&"bob", &"carol", 2)?;
|
|
//! history.record_winner(&"alice", &"carol", 3)?;
|
|
//!
|
|
//! let report = history.converge()?;
|
|
//! assert!(report.converged);
|
|
//!
|
|
//! let alice = history.current_skill("alice").unwrap();
|
|
//! assert!(alice.mu() > 0.0, "alice won every game she played");
|
|
//! # Ok::<(), trueskill_tt::InferenceError>(())
|
|
//! ```
|
|
//!
|
|
//! Teams, weights, explicit rankings and continuous scores go through the
|
|
//! fluent event builder:
|
|
//!
|
|
//! ```
|
|
//! use trueskill_tt::History;
|
|
//!
|
|
//! let mut history = History::builder().p_draw(0.1).build();
|
|
//!
|
|
//! history
|
|
//! .event(1)
|
|
//! .team(["alice", "bob"])
|
|
//! .team(["carol", "dave"])
|
|
//! .ranking([0, 1])
|
|
//! .commit()?;
|
|
//!
|
|
//! history.converge()?;
|
|
//! # Ok::<(), trueskill_tt::InferenceError>(())
|
|
//! ```
|
|
//!
|
|
//! # Draws need a draw probability
|
|
//!
|
|
//! A `p_draw` of zero asserts that draws cannot happen, so a tied result has
|
|
//! no representable likelihood and is rejected:
|
|
//!
|
|
//! ```
|
|
//! use trueskill_tt::{History, InferenceError};
|
|
//!
|
|
//! let mut history = History::default(); // p_draw defaults to 0.0
|
|
//! let err = history.record_draw(&"alice", &"bob", 1).unwrap_err();
|
|
//! assert!(matches!(err, InferenceError::TieWithoutDrawProbability { .. }));
|
|
//! ```
|
|
//!
|
|
//! This also applies to [`Outcome::winner`] for three or more teams, which
|
|
//! ties every loser. Configure a positive `p_draw` for those.
|
|
//!
|
|
//! # Core types
|
|
//!
|
|
//! - [`History`] — the top-level container: ingests events, runs
|
|
//! forward/backward message passing, and answers queries.
|
|
//! - [`Gaussian`] — the probability type, stored in natural parameters
|
|
//! (`pi = 1/sigma²`, `tau = mu/sigma²`) so message passing is add/subtract.
|
|
//! - [`Game`] — one match in isolation, for scoring a hypothetical without a
|
|
//! history.
|
|
//! - [`Outcome`] — how a match ended: ranks, or continuous scores.
|
|
//! - [`Rating`] — a competitor's static configuration (prior, `beta`, drift).
|
|
//!
|
|
//! # Feature flags
|
|
//!
|
|
//! - `approx` — implements [`approx`](https://docs.rs/approx) equality traits
|
|
//! for [`Gaussian`]. Useful in tests.
|
|
//! - `rayon` — parallelises the within-slice sweep and the per-slice passes of
|
|
//! `learning_curves`/`log_evidence`. Opt-in; results stay bit-identical
|
|
//! regardless of worker count.
|
|
|
|
#![forbid(unsafe_code)]
|
|
// Turned on once the surface was fully documented (80 items at the time), so
|
|
// the next undocumented public item is a build failure rather than a warning
|
|
// nobody reads.
|
|
#![deny(missing_docs)]
|
|
|
|
/// Compiles every `rust` block in `README.md` as a doctest.
|
|
///
|
|
/// The README is not the crate's front page — the module docs above are — so it
|
|
/// is pulled in here rather than via a crate-level `#![doc = ...]`, purely so
|
|
/// its examples are type-checked. Without this nothing compiled them, and they
|
|
/// had drifted far enough that four blocks no longer built (#35). `cfg(doctest)`
|
|
/// means this type exists only while collecting doctests.
|
|
///
|
|
/// Blocks that are illustrative rather than runnable are fenced as `text`.
|
|
#[cfg(doctest)]
|
|
#[doc = include_str!("../README.md")]
|
|
pub struct ReadmeDoctests;
|
|
|
|
use std::{
|
|
cmp::Reverse,
|
|
f64::consts::{FRAC_1_SQRT_2, FRAC_2_SQRT_PI, SQRT_2},
|
|
};
|
|
|
|
mod acquisition;
|
|
#[cfg(feature = "approx")]
|
|
mod approx;
|
|
pub(crate) mod arena;
|
|
mod color_group;
|
|
mod competitor;
|
|
mod convergence;
|
|
/// Skill drift: how much a competitor's skill is allowed to move between
|
|
/// appearances.
|
|
///
|
|
/// Public because [`Drift`] is a trait a caller may implement — a per-sport
|
|
/// off-season, say, or a schedule where drift is a function of the calendar
|
|
/// rather than of elapsed ticks. [`ConstantDrift`] is what
|
|
/// [`HistoryBuilder`] uses by default.
|
|
pub mod drift;
|
|
mod error;
|
|
mod event;
|
|
mod event_builder;
|
|
pub(crate) mod factor;
|
|
mod game;
|
|
/// The Gaussian message type and its expectation-propagation algebra.
|
|
///
|
|
/// Public because [`Gaussian`] appears throughout the results: a posterior
|
|
/// skill, a learning-curve point, a predicted margin. The module carries the
|
|
/// operator documentation — `Mul`/`Div` are the EP product and cavity, not
|
|
/// arithmetic on random variables.
|
|
pub mod gaussian;
|
|
mod history;
|
|
mod joint;
|
|
mod key_table;
|
|
mod matrix;
|
|
mod observer;
|
|
mod outcome;
|
|
mod predict;
|
|
pub(crate) mod quadrature;
|
|
mod rating;
|
|
pub(crate) mod storage;
|
|
mod time;
|
|
mod time_slice;
|
|
|
|
pub use acquisition::expected_information_gain;
|
|
pub use convergence::{ConvergenceOptions, ConvergenceReport};
|
|
pub use drift::{ConstantDrift, Drift};
|
|
pub use error::{InferenceError, UnknownKeys};
|
|
pub use event::{Event, Member, Team};
|
|
pub use event_builder::EventBuilder;
|
|
pub use game::{Game, GameOptions};
|
|
pub use gaussian::Gaussian;
|
|
pub use history::{History, HistoryBuilder, Joint};
|
|
use matrix::Matrix;
|
|
pub use observer::{NullObserver, Observer};
|
|
pub use outcome::Outcome;
|
|
pub use predict::Prediction;
|
|
pub use rating::Rating;
|
|
/// The `smallvec` crate, re-exported.
|
|
///
|
|
/// Four public items name `SmallVec` in their signatures: [`Event::teams`],
|
|
/// [`Team::members`], [`Outcome::Ranked`]'s payload and
|
|
/// [`ConvergenceReport::per_iteration_time`]. You can *build* an `Event`
|
|
/// without ever naming the type — `vec![..].into()` and `.collect()` both work
|
|
/// — and iterate the timings through `Deref`. But writing a helper that
|
|
/// *returns* a teams list, or a `match` arm that binds ranks and passes them
|
|
/// on, requires the type by name.
|
|
///
|
|
/// Measured: the only `Joint` doc example failed to compile from a consumer
|
|
/// crate with `unresolved import \`smallvec\``, because the dependency was in
|
|
/// the signature but not reachable. Re-exported so a consumer takes this
|
|
/// crate's version rather than pinning a matching one of their own.
|
|
pub use smallvec;
|
|
pub use time::{Time, Untimed};
|
|
|
|
/// Default performance noise: how much a single showing varies around skill.
|
|
///
|
|
/// Every other default is expressed as a multiple of this, so `BETA` sets the
|
|
/// scale of the whole rating system. Doubling it and doubling `SIGMA` and
|
|
/// `GAMMA` with it gives the same fit on a rescaled axis.
|
|
pub const BETA: f64 = 1.0;
|
|
/// Default prior mean skill.
|
|
///
|
|
/// Zero rather than a conventional 25: the scale is set by `BETA`, and a
|
|
/// centred axis makes a negative rating mean "below the prior" instead of
|
|
/// looking like an error.
|
|
pub const MU: f64 = 0.0;
|
|
/// Default prior standard deviation: how unsure the model starts out.
|
|
///
|
|
/// Six betas is deliberately wide — a new competitor's first result should
|
|
/// move them a long way, and the prior should not fight the evidence.
|
|
pub const SIGMA: f64 = BETA * 6.0;
|
|
/// Default drift: the standard deviation of skill movement per unit of time.
|
|
///
|
|
/// Enters inference as a *variance* (`gamma^2` per elapsed tick), which is why
|
|
/// [`ConstantDrift`] squares it and why a negative gamma would be
|
|
/// indistinguishable from its absolute value — see [`ConstantDrift::new`].
|
|
pub const GAMMA: f64 = BETA * 0.03;
|
|
/// Default draw probability: zero, meaning ties are not modelled.
|
|
///
|
|
/// A history that ingests a tie needs a positive value. With `p_draw == 0.0`
|
|
/// the truncation margin is zero and the two-sided tie update evaluates
|
|
/// `0/0`, so ingestion rejects such events with
|
|
/// [`InferenceError::TieWithoutDrawProbability`].
|
|
pub const P_DRAW: f64 = 0.0;
|
|
/// Default convergence threshold, in the same units as
|
|
/// [`ConvergenceReport::final_step`](crate::ConvergenceReport).
|
|
///
|
|
/// The sweep stops once the largest change a full iteration makes to any
|
|
/// message falls below this.
|
|
pub const EPSILON: f64 = 1e-6;
|
|
/// Default cap on convergence sweeps.
|
|
///
|
|
/// **A runaway guard, not a budget.** The sweep exits as soon as the step falls
|
|
/// below `epsilon`, so the cap is never reached by a history that converges and
|
|
/// raising it costs nothing. Measured on a history that needs four sweeps:
|
|
///
|
|
/// ```text
|
|
/// max_iter 30: 4 iterations, 129.9 us
|
|
/// max_iter 100_000: 4 iterations, 131.9 us
|
|
/// ```
|
|
///
|
|
/// This was `30` until it was measured, and 30 truncated ordinary healthy
|
|
/// histories: 160 events over 100 competitors already needs 42. Because a short
|
|
/// fit is finite and sensibly ordered, that was invisible.
|
|
///
|
|
/// # Why it is not scaled to the history
|
|
///
|
|
/// The obvious improvement — pick the cap from the node or event count — does
|
|
/// not work, because iteration count is driven by how *loopy* the graph is
|
|
/// rather than how big it is. At a fixed 320 events over 40 slices, varying
|
|
/// only the number of competitors sharing them:
|
|
///
|
|
/// ```text
|
|
/// competitors appearances each iterations
|
|
/// 3 213 2_789
|
|
/// 10 64 1_068
|
|
/// 50 12.8 206
|
|
/// 100 6.4 90
|
|
/// 400 1.6 2
|
|
/// ```
|
|
///
|
|
/// Three orders of magnitude apart on identical event and slice counts. Any
|
|
/// formula in those two numbers would be badly wrong on some real shape, so the
|
|
/// cap is a single value set high enough that reaching it means the fit is
|
|
/// oscillating rather than merely large.
|
|
///
|
|
/// Reaching it is [`InferenceError::NotConverged`]. See
|
|
/// [`History::converge`](crate::History::converge).
|
|
pub const ITERATIONS: usize = 10_000;
|
|
|
|
/// 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;
|
|
/// `1 / sqrt(pi)`, the leading factor of the `erfcx` continued fraction.
|
|
const FRAC_1_SQRT_PI: f64 = 0.564_189_583_547_756_3;
|
|
/// `sqrt(2 / pi)`, the numerator of the inverse Mills ratio in scaled form.
|
|
const SQRT_2_OVER_PI: f64 = 0.797_884_560_802_865_4;
|
|
/// How many window widths into the tail before a tie window is treated as a
|
|
/// half-line. Beyond this the truncated mass is concentrated within `1/alpha`
|
|
/// of the near edge, so the far edge contributes nothing measurable.
|
|
const HALF_LINE_WINDOW: f64 = 10.0;
|
|
/// Where `v - alpha` switches from subtraction to its asymptotic series.
|
|
///
|
|
/// The subtraction loses roughly `eps * alpha^2` of relative precision, and the
|
|
/// four-term series is good to ~1e-10 by here, so the two are at their closest
|
|
/// agreement around this point. Below it the subtraction is exact; above it the
|
|
/// series is.
|
|
/// `alpha / width` past which the tie branch's `v^2 - u` has lost too many
|
|
/// digits to trust, and the narrow-window form takes over.
|
|
///
|
|
/// The subtraction retains about `(width / alpha)^2 / EPSILON` of its
|
|
/// precision, so this is the ratio at which that falls below roughly 1e-6.
|
|
const NARROW_WINDOW_RATIO: f64 = 2.0e4;
|
|
const ASYMPTOTIC_MILLS_ALPHA: f64 = 100.0;
|
|
|
|
pub(crate) const N00: Gaussian = Gaussian::from_ms(0.0, 0.0);
|
|
pub(crate) const N_INF: Gaussian = Gaussian::from_ms(0.0, f64::INFINITY);
|
|
|
|
/// An interned competitor handle: a dense slot number, not a user key.
|
|
///
|
|
/// [`History`] stores skills and messages by `Index` rather than by `K`, so
|
|
/// the hot path never hashes a key. [`History::intern`] promotes a key to one
|
|
/// and [`History::lookup`] resolves an existing key without creating.
|
|
///
|
|
/// Indices are assigned in interning order and are stable for the life of a
|
|
/// history. They are **not** portable between histories: the same key interns
|
|
/// to different slots depending on ingestion order.
|
|
#[derive(Copy, Clone, Default, PartialEq, PartialOrd, Eq, Ord, Hash, Debug)]
|
|
pub struct Index(usize);
|
|
|
|
impl Index {
|
|
/// The underlying slot number.
|
|
///
|
|
/// Indices are dense and assigned in interning order, so this is usable as
|
|
/// a key into a caller-side side table.
|
|
#[must_use]
|
|
pub fn get(self) -> usize {
|
|
self.0
|
|
}
|
|
}
|
|
|
|
impl From<usize> for Index {
|
|
fn from(ix: usize) -> Self {
|
|
Self(ix)
|
|
}
|
|
}
|
|
|
|
impl From<Index> for usize {
|
|
fn from(idx: Index) -> Self {
|
|
idx.0
|
|
}
|
|
}
|
|
|
|
/// Complementary error function.
|
|
///
|
|
/// # Why every transcendental in this crate goes through `libm`
|
|
///
|
|
/// IEEE 754 specifies the basic operations and `sqrt` exactly, but says nothing
|
|
/// about `exp`, `log` or `erf`. `std`'s versions delegate to the *system* math
|
|
/// library, so they differ between platforms: measured here, `f64::exp` and
|
|
/// `libm::exp` disagree on 9.7% of inputs and `f64::ln` / `libm::log` on 5.0%,
|
|
/// each by one ULP.
|
|
///
|
|
/// Inference is an iterative fixed point, so a one-ULP difference can change an
|
|
/// iteration count and therefore the answer by more than one ULP. Routing every
|
|
/// transcendental through `libm` makes a fit reproducible across platforms, not
|
|
/// just across thread counts as `tests/determinism.rs` already checks.
|
|
///
|
|
/// **So: use `libm::exp` / `libm::log` in inference code, never `f64::exp` /
|
|
/// `f64::ln`.** `sqrt` is exempt — IEEE specifies it exactly, so `f64::sqrt` is
|
|
/// already portable. Test code may use whichever is clearer.
|
|
///
|
|
/// It costs nothing: `Batch::iteration` measured -2.7% [-5.7%, -0.3%] with the
|
|
/// whole set swapped.
|
|
///
|
|
/// Delegates to `libm`, which is the Rust port of FDLIBM and accurate to about
|
|
/// one ULP. This replaced a Numerical Recipes `erfcc` rational approximation
|
|
/// whose documented bound was 1.2e-7 *relative* — measured at ~1e-7 across the
|
|
/// whole range, and the binding accuracy constraint on the entire crate.
|
|
///
|
|
/// The swap is free. 98% of the arguments inference passes here have
|
|
/// `|x| < 0.84375`, which is exactly where FDLIBM skips the exponential
|
|
/// entirely, so the longer polynomial costs nothing on the distribution that
|
|
/// actually occurs: `Batch::iteration` moved -1.6% [-4.7%, +0.9%], p = 0.31.
|
|
///
|
|
/// What it bought: `compute_margin` went from 8.4e-8 to 1.7e-16 against exact
|
|
/// quantiles, `cdf(mu, mu, sigma)` is now exactly 0.5, and `sf + cdf` sums to
|
|
/// one within a single ULP where it was 3e-8 out.
|
|
fn erfc(x: f64) -> f64 {
|
|
libm::erfc(x)
|
|
}
|
|
|
|
/// The previous Numerical Recipes `erfcc`, kept only so the timing test can
|
|
/// compare both in one binary. Removed once the comparison is recorded.
|
|
fn erfc_inv(mut y: f64) -> f64 {
|
|
if y >= 2.0 {
|
|
return f64::NEG_INFINITY;
|
|
}
|
|
|
|
debug_assert!(y >= 0.0, "y must be nonnegative");
|
|
|
|
if y == 0.0 {
|
|
return f64::INFINITY;
|
|
}
|
|
|
|
if y >= 1.0 {
|
|
y = 2.0 - y;
|
|
}
|
|
|
|
let t = libm::sqrt(-2.0 * libm::log(y / 2.0));
|
|
|
|
// The leading coefficient is NEGATIVE. `rational - t` is negative here, so
|
|
// a positive coefficient mirrors the starting point to `-x0` — the
|
|
// reflection of the root. Newton then has to cross the origin to get back,
|
|
// which a fixed iteration count does not manage: measured against the true
|
|
// value, `erfc_inv(0.1)` returned 1.044 instead of 1.16309, and the error
|
|
// grew as y shrank until `compute_margin` stopped being monotone in
|
|
// `p_draw` altogether.
|
|
let mut x =
|
|
-FRAC_1_SQRT_2 * ((2.30753 + t * 0.27061) / (1.0 + t * (0.99229 + t * 0.04481)) - t);
|
|
|
|
for _ in 0..3 {
|
|
let err = erfc(x) - y;
|
|
|
|
x += err / (FRAC_2_SQRT_PI * libm::exp(-(x * x)) - x * err)
|
|
}
|
|
|
|
if y < 1.0 { x } else { -x }
|
|
}
|
|
|
|
fn ppf(p: f64, mu: f64, sigma: f64) -> f64 {
|
|
mu - sigma * SQRT_2 * erfc_inv(2.0 * p)
|
|
}
|
|
|
|
fn compute_margin(p_draw: f64, sd: f64) -> f64 {
|
|
ppf(0.5 - p_draw / 2.0, 0.0, sd).abs()
|
|
}
|
|
|
|
pub(crate) fn cdf(x: f64, mu: f64, sigma: f64) -> f64 {
|
|
let z = -(x - mu) / (sigma * SQRT_2);
|
|
|
|
0.5 * erfc(z)
|
|
}
|
|
|
|
/// `P(X > x)` for `X ~ N(mu, sigma^2)`.
|
|
///
|
|
/// The survival function, computed directly rather than as `1 - cdf(..)`.
|
|
///
|
|
/// The two are algebraically identical and numerically are not. `cdf` returns
|
|
/// a value approaching 1 for an upper tail, so subtracting it from 1 cancels
|
|
/// away every significant digit the tail had: measured against this function,
|
|
/// `1 - cdf` carries 7% error by four sigma past the mean and returns exactly
|
|
/// zero beyond about 8.3 sigma — where the true value is still 1e-19 and
|
|
/// perfectly representable. `erfc` holds *relative* accuracy all the way down
|
|
/// to 1e-296, so the precision is there to keep; only the subtraction threw it
|
|
/// away.
|
|
///
|
|
/// This matters most where evidence is smallest, which is exactly where an
|
|
/// upset makes it interesting: `ln` of a clamped zero is -708 regardless of
|
|
/// whether the truth was -43 or -600.
|
|
pub(crate) fn sf(x: f64, mu: f64, sigma: f64) -> f64 {
|
|
0.5 * erfc((x - mu) / (sigma * SQRT_2))
|
|
}
|
|
|
|
/// `e^(x^2) * erfc(x)`, the scaled complementary error function, for `x >= 0`.
|
|
///
|
|
/// Exists so the exponential factor common to a Gaussian density and its tail
|
|
/// integral can be cancelled *analytically* instead of being computed twice
|
|
/// and divided. Both underflow to zero past about 26 sigma, and their ratio is
|
|
/// then `0/0` — finite in the limit, `NaN` in floating point.
|
|
fn erfcx(x: f64) -> f64 {
|
|
if x < 2.0 {
|
|
// Below the crossover neither factor is extreme: erfc is O(1) and
|
|
// exp(x^2) is at most e^4, so the direct product is exact enough and
|
|
// cheaper than the continued fraction.
|
|
libm::exp(x * x) * erfc(x)
|
|
} else {
|
|
// erfcx(x) = 1/sqrt(pi) * 1/(x + (1/2)/(x + 1/(x + (3/2)/(x + ...)))),
|
|
// evaluated by backward recurrence. Converges quickly for x >= 2 and,
|
|
// unlike the product form, never touches an exponential.
|
|
let mut f = 0.0;
|
|
for n in (1..=60u32).rev() {
|
|
f = (f64::from(n) * 0.5) / (x + f);
|
|
}
|
|
FRAC_1_SQRT_PI / (x + f)
|
|
}
|
|
}
|
|
|
|
/// `ln` of the normal density at `x`.
|
|
///
|
|
/// The density itself underflows to zero past about 38 sigma, and `ln` of a
|
|
/// clamped zero is -708 whatever the truth was. The log form is a polynomial:
|
|
/// it stays exact at any separation, and the values it produces (-5001 nats at
|
|
/// 100 sigma, -500001 at 1000) are perfectly representable.
|
|
pub(crate) fn ln_pdf(x: f64, mu: f64, sigma: f64) -> f64 {
|
|
let z = (x - mu) / sigma;
|
|
-libm::log(SQRT_TAU * sigma) - 0.5 * z * z
|
|
}
|
|
|
|
/// `ln P(X > x)` for `X ~ N(mu, sigma^2)`.
|
|
///
|
|
/// In the upper tail the `exp(-z^2 / 2)` common to the tail integral is
|
|
/// factored out analytically via `erfcx`, so this never underflows — where
|
|
/// `sf(..).ln()` bottoms out at -708 once `erfc` itself reaches zero.
|
|
pub(crate) fn ln_sf(x: f64, mu: f64, sigma: f64) -> f64 {
|
|
let z = (x - mu) / sigma;
|
|
|
|
if z > 0.0 {
|
|
// ln(0.5 * erfc(z/sqrt2)) with erfc(y) = exp(-y^2) * erfcx(y).
|
|
-std::f64::consts::LN_2 - 0.5 * z * z + libm::log(erfcx(z / SQRT_2))
|
|
} else {
|
|
// The mass here is at least a half; nothing to lose.
|
|
libm::log(sf(x, mu, sigma))
|
|
}
|
|
}
|
|
|
|
/// `ln P(lo < X < hi)` for `X ~ N(mu, sigma^2)`.
|
|
///
|
|
/// When the interval sits in a tail both endpoint probabilities underflow
|
|
/// together, so their difference is taken in scaled form with the shared
|
|
/// exponential factored out. When it straddles the mean nothing is small and
|
|
/// the direct difference is exact.
|
|
pub(crate) fn ln_interval(lo: f64, hi: f64, mu: f64, sigma: f64) -> f64 {
|
|
let z_lo = (lo - mu) / sigma;
|
|
let z_hi = (hi - mu) / sigma;
|
|
|
|
if z_hi <= z_lo {
|
|
return f64::NEG_INFINITY;
|
|
}
|
|
|
|
// Fold a lower-tail interval onto the upper tail; the normal is symmetric.
|
|
let (near, far) = if z_lo >= 0.0 {
|
|
(z_lo, z_hi)
|
|
} else if z_hi <= 0.0 {
|
|
(-z_hi, -z_lo)
|
|
} else {
|
|
// Straddles the mean: the interval holds a non-negligible share of the
|
|
// mass, so neither endpoint is near enough to 1 to cancel.
|
|
return libm::log((cdf(hi, mu, sigma) - cdf(lo, mu, sigma)).max(f64::MIN_POSITIVE));
|
|
};
|
|
|
|
let (a, b) = (near / SQRT_2, far / SQRT_2);
|
|
// b > a >= 0, so this ratio of exponentials is at most 1 and cannot overflow.
|
|
let scale = libm::exp(a * a - b * b);
|
|
let bracket = erfcx(a) - scale * erfcx(b);
|
|
|
|
if bracket <= 0.0 {
|
|
return f64::NEG_INFINITY;
|
|
}
|
|
|
|
-std::f64::consts::LN_2 - a * a + libm::log(bracket)
|
|
}
|
|
|
|
fn pdf(x: f64, mu: f64, sigma: f64) -> f64 {
|
|
let normalizer = (SQRT_TAU * sigma).powi(-1);
|
|
let functional = libm::exp(-((x - mu) * (x - mu)) / (2.0 * sigma * sigma));
|
|
|
|
normalizer * functional
|
|
}
|
|
|
|
/// Truncated-Gaussian correction terms `(v, w)`.
|
|
///
|
|
/// `v` shifts the mean and `w` shrinks the variance. Both are ratios whose
|
|
/// numerator and denominator underflow together in the tails, so both are
|
|
/// computed in scaled form there: the shared `exp(-alpha^2 / 2)` is cancelled
|
|
/// analytically rather than evaluated and divided out. Without that, a
|
|
/// truncation point beyond about 39 sigma produced `0 / 0` and put `NaN`
|
|
/// straight into the posterior.
|
|
/// Truncation terms for a boundary `alpha` standard deviations into the upper
|
|
/// tail, from the asymptotic expansion of the inverse Mills ratio.
|
|
///
|
|
/// `v` tends to `alpha` out here, so the gap between them cannot be obtained by
|
|
/// subtracting one from the other — the series computes the gap directly, and
|
|
/// `w = v * gap` then never forms the difference of two large near-equal
|
|
/// numbers. A far-tail *window* behaves like a half-line once it is more than a
|
|
/// few multiples of its own width from the mean, so the tie branch shares this.
|
|
fn half_line_truncation(alpha: f64) -> (f64, f64) {
|
|
let inv = alpha.recip();
|
|
let inv_sq = inv * inv;
|
|
let b = 2.0 - inv_sq * (10.0 - 74.0 * inv_sq);
|
|
let gap = inv * (1.0 - inv_sq * b);
|
|
let v = alpha + gap;
|
|
|
|
// Returns `1 - w`, not `w`, and that is the whole point of this shape.
|
|
//
|
|
// `w` tends to 1 out here, so a caller forming `1 - w` loses about
|
|
// `log10(alpha^2)` digits: measured against the exact truncated variance,
|
|
// `1 - w` came back with 8.9e-5 relative error at alpha = 1e6 and **0.0**
|
|
// from alpha = 1e8 — where the true value is 1e-16 and perfectly
|
|
// representable. `sigma * (1 - w).sqrt()` was then exactly zero, and
|
|
// `from_ms(mu, 0.0)` is a point mass whose `mu()` is `inf/inf = NaN`.
|
|
//
|
|
// Expanding `1 - v*gap` symbolically removes the subtraction: with
|
|
// `alpha*gap = 1 - inv^2*b`, the leading ones cancel on paper instead of in
|
|
// floating point, leaving `inv^2` times a bracket that tends to 1. Measured
|
|
// exact — 0.0 relative error — from alpha = 1e3 to 1e8.
|
|
let one_minus_w = inv_sq
|
|
* ((1.0 - inv_sq * (10.0 - 74.0 * inv_sq)) + 2.0 * inv_sq * b - inv_sq * inv_sq * b * b);
|
|
|
|
(v, one_minus_w)
|
|
}
|
|
|
|
/// Truncation to a *narrow* window `[alpha, alpha + d]`, as `(v, 1 - w)`.
|
|
///
|
|
/// The tie branch forms `w` from `v^2 - u`, and both grow as `alpha^2` while
|
|
/// their difference stays `O(1)`. Far enough into the tail that subtraction has
|
|
/// nothing left: measured at `alpha = 1e6` with a window of `1e-6` it kept four
|
|
/// significant digits and returned `1 - w = -2.4e-4` where the truth is
|
|
/// `+2.8e-13`, so `sqrt` of it was NaN. One step earlier it was quietly wrong
|
|
/// instead — `1 - w = 1.0` exactly, a truncation reported as a no-op, where the
|
|
/// truth was `5e-17`.
|
|
///
|
|
/// The existing half-line escape hatch does not cover it, because that keys on
|
|
/// `alpha * d >= HALF_LINE_WINDOW` — how many window-widths from the mean the
|
|
/// window sits — and a *narrow* window fails that however deep it is.
|
|
///
|
|
/// Over a narrow window the density is `exp(-t*s - s^2 d^2 / 2)` in
|
|
/// `x = alpha + s*d`, with `t = alpha * d`. Dropping the `d^2` term leaves a
|
|
/// truncated exponential on `[0, 1]`, whose mean and variance are closed forms.
|
|
/// So `v = alpha + d*m(t)` and `1 - w = d^2 * V(t)`, with no subtraction of
|
|
/// large quantities anywhere.
|
|
///
|
|
/// Measured against high-precision quadrature over `alpha` in `[1e2, 1e9]`:
|
|
/// `v` exact to 4e-10 or better, `1 - w` to 4e-10 across the region this is
|
|
/// used in.
|
|
fn narrow_window_truncation(alpha: f64, d: f64) -> (f64, f64) {
|
|
let t = alpha * d;
|
|
|
|
// `m` and `V` are the mean and variance of a truncated exponential on
|
|
// [0, 1] with rate `t`, both of which cancel as `t -> 0`. The series is
|
|
// their limit (1/2 and 1/12, a uniform window) with the leading correction.
|
|
let (m, v_s) = if t < 1e-3 {
|
|
(
|
|
0.5 - t / 12.0 + t * t * t / 720.0,
|
|
1.0 / 12.0 - t * t / 240.0,
|
|
)
|
|
} else {
|
|
let em1 = libm::expm1(t);
|
|
(
|
|
1.0 / t - 1.0 / em1,
|
|
1.0 / (t * t) - (em1 + 1.0) / (em1 * em1),
|
|
)
|
|
};
|
|
|
|
(alpha + d * m, d * d * v_s)
|
|
}
|
|
|
|
fn v_w(mu: f64, sigma: f64, margin: f64, tie: bool) -> (f64, f64) {
|
|
if !tie {
|
|
let alpha = (margin - mu) / sigma;
|
|
|
|
// v is the inverse Mills ratio, phi(alpha) / Phi(-alpha), and w needs
|
|
// the gap `v - alpha` as well as v itself. Far into the tail v tends to
|
|
// alpha, so that gap is a subtraction of two nearly equal numbers and
|
|
// loses every digit it has: at alpha = 1e6 it drove w above 1 and made
|
|
// `sqrt(1 - w)` NaN. Past the crossover the gap comes from its
|
|
// asymptotic series instead, which has no subtraction in it.
|
|
if alpha >= ASYMPTOTIC_MILLS_ALPHA {
|
|
return half_line_truncation(alpha);
|
|
}
|
|
|
|
let (v, gap) = if alpha > 0.0 {
|
|
// Both terms carry exp(-alpha^2 / 2); in scaled form it cancels
|
|
// and the result stays exact however far into the tail alpha sits.
|
|
let v = SQRT_2_OVER_PI / erfcx(alpha / SQRT_2);
|
|
(v, v - alpha)
|
|
} else {
|
|
// Phi(-alpha) >= 1/2 here, so the direct ratio loses nothing.
|
|
let v = pdf(-alpha, 0.0, 1.0) / cdf(-alpha, 0.0, 1.0);
|
|
(v, v - alpha)
|
|
};
|
|
|
|
(v, 1.0 - v * gap)
|
|
} else {
|
|
// v is odd in mu and w is even, so fold to mu <= 0. Both truncation
|
|
// points then sit in the upper tail, where the scaled form applies.
|
|
let flipped = mu > 0.0;
|
|
let mu = if flipped { -mu } else { mu };
|
|
|
|
let alpha = (-margin - mu) / sigma;
|
|
let beta = (margin - mu) / sigma;
|
|
|
|
// `w` comes out of `v * v - u`, and both terms grow as alpha^2 while
|
|
// their difference stays O(1) — at alpha = 1e9 that subtraction had no
|
|
// digits left and returned w = -128, making `sqrt(1 - w)` nonsense.
|
|
// Once the window sits many of its own widths into the tail it is
|
|
// indistinguishable from a half-line, so the asymptotic covers it with
|
|
// no subtraction at all.
|
|
let width = beta - alpha;
|
|
|
|
if alpha >= ASYMPTOTIC_MILLS_ALPHA && alpha * width >= HALF_LINE_WINDOW {
|
|
let (v, one_minus_w) = half_line_truncation(alpha);
|
|
return (if flipped { -v } else { v }, one_minus_w);
|
|
}
|
|
|
|
// A narrow window deep in the tail: too narrow for the half-line above,
|
|
// too deep for the subtraction below. The direct form keeps roughly
|
|
// `1 / (alpha/width)^2` of its digits, so the crossover is on that
|
|
// ratio rather than on either quantity alone — and the approximation is
|
|
// most accurate exactly where the subtraction is worst, since both
|
|
// improve as the window narrows.
|
|
if alpha > 0.0 && alpha > NARROW_WINDOW_RATIO * width {
|
|
let (v, one_minus_w) = narrow_window_truncation(alpha, width);
|
|
return (if flipped { -v } else { v }, one_minus_w);
|
|
}
|
|
|
|
let (v, u) = if alpha > 0.0 {
|
|
// beta > alpha > 0, so this ratio of exponentials is at most 1 and
|
|
// cannot overflow.
|
|
let scale = libm::exp(0.5 * (alpha * alpha - beta * beta));
|
|
let denominator = 0.5 * (erfcx(alpha / SQRT_2) - scale * erfcx(beta / SQRT_2));
|
|
|
|
(
|
|
(1.0 - scale) / SQRT_TAU / denominator,
|
|
(alpha - beta * scale) / SQRT_TAU / denominator,
|
|
)
|
|
} else {
|
|
// The interval straddles the mean, so nothing here is small.
|
|
let denominator = cdf(beta, 0.0, 1.0) - cdf(alpha, 0.0, 1.0);
|
|
|
|
(
|
|
(pdf(alpha, 0.0, 1.0) - pdf(beta, 0.0, 1.0)) / denominator,
|
|
(alpha * pdf(alpha, 0.0, 1.0) - beta * pdf(beta, 0.0, 1.0)) / denominator,
|
|
)
|
|
};
|
|
|
|
// `1 - w` where `w = v^2 - u`. Both `v^2` and `u` grow as alpha^2 while
|
|
// their difference stays O(1), so this subtraction is the one place the
|
|
// tie branch can still lose everything — see the escape hatch above,
|
|
// which is what keeps the far tail away from it.
|
|
let one_minus_w = 1.0 + u - v.powi(2);
|
|
|
|
(if flipped { -v } else { v }, one_minus_w)
|
|
}
|
|
}
|
|
|
|
fn trunc(mu: f64, sigma: f64, margin: f64, tie: bool) -> (f64, f64) {
|
|
// `v_w` returns `1 - w` rather than `w`: forming the difference here is
|
|
// what destroyed the truncated variance in the far tail.
|
|
let (v, one_minus_w) = v_w(mu, sigma, margin, tie);
|
|
|
|
let mu_trunc = mu + sigma * v;
|
|
let sigma_trunc = sigma * one_minus_w.sqrt();
|
|
|
|
(mu_trunc, sigma_trunc)
|
|
}
|
|
|
|
pub(crate) fn approx(n: Gaussian, margin: f64, tie: bool) -> Gaussian {
|
|
let (mu, sigma) = trunc(n.mu(), n.sigma(), margin, tie);
|
|
|
|
Gaussian::from_ms(mu, sigma)
|
|
}
|
|
|
|
/// Componentwise maximum that **propagates** NaN rather than dropping it.
|
|
///
|
|
/// Every caller folds this as `tuple_max(accumulator, new)`. A plain `>`
|
|
/// comparison is false against NaN, so a NaN accumulator would be replaced by
|
|
/// the next finite delta and the breakdown would vanish — leaving `step_is_finite`
|
|
/// to pass on a fit that is already NaN. Because the fold runs over a `HashMap`,
|
|
/// whether that happened depended on per-process hash order: measured, a NaN fit
|
|
/// was reported as `converged: true` in 16 of 30 runs on identical input.
|
|
///
|
|
/// `f64::max` is not a substitute: it also ignores NaN by design, which is the
|
|
/// same defect wearing a standard-library name.
|
|
pub(crate) fn tuple_max(v1: (f64, f64), v2: (f64, f64)) -> (f64, f64) {
|
|
(
|
|
max_propagating_nan(v1.0, v2.0),
|
|
max_propagating_nan(v1.1, v2.1),
|
|
)
|
|
}
|
|
|
|
fn max_propagating_nan(a: f64, b: f64) -> f64 {
|
|
if a.is_nan() || b.is_nan() {
|
|
f64::NAN
|
|
} else if a > b {
|
|
a
|
|
} else {
|
|
b
|
|
}
|
|
}
|
|
|
|
pub(crate) fn tuple_gt(t: (f64, f64), e: f64) -> bool {
|
|
t.0 > e || t.1 > e
|
|
}
|
|
|
|
/// Whether a convergence step is finite in both components.
|
|
///
|
|
/// A NaN step means EP broke down numerically. Because every comparison
|
|
/// against NaN is false, `tuple_gt` reads NaN as "below epsilon" — so
|
|
/// convergence checks must test finiteness explicitly rather than inferring
|
|
/// success from `!tuple_gt(..)`.
|
|
pub(crate) fn step_is_finite(t: (f64, f64)) -> bool {
|
|
t.0.is_finite() && t.1.is_finite()
|
|
}
|
|
|
|
/// Whether a step counts as converged: finite *and* within `epsilon`.
|
|
pub(crate) fn step_converged(t: (f64, f64), epsilon: f64) -> bool {
|
|
step_is_finite(t) && !tuple_gt(t, epsilon)
|
|
}
|
|
|
|
/// Indices of the first pair of teams sharing a rank, if any.
|
|
///
|
|
/// A tie is only representable when the draw probability is positive: with
|
|
/// `p_draw == 0.0` the truncation margin collapses to zero and the two-sided
|
|
/// tie update evaluates `0/0`. Callers use this to reject such events before
|
|
/// they reach inference.
|
|
pub(crate) fn first_tied_pair(ranks: &[u32]) -> Option<(usize, usize)> {
|
|
for (i, a) in ranks.iter().enumerate() {
|
|
for (j, b) in ranks.iter().enumerate().skip(i + 1) {
|
|
if a == b {
|
|
return Some((i, j));
|
|
}
|
|
}
|
|
}
|
|
|
|
None
|
|
}
|
|
|
|
/// As `first_tied_pair`, but over the engine's internal `f64` outputs.
|
|
///
|
|
/// Ranks reach the engine already converted to descending `f64` outputs, and
|
|
/// `Game` decides a tie by exact equality of those values — so this mirrors
|
|
/// the comparison inference itself performs.
|
|
pub(crate) fn first_tied_output(outputs: &[f64]) -> Option<(usize, usize)> {
|
|
for (i, a) in outputs.iter().enumerate() {
|
|
for (j, b) in outputs.iter().enumerate().skip(i + 1) {
|
|
if a == b {
|
|
return Some((i, j));
|
|
}
|
|
}
|
|
}
|
|
|
|
None
|
|
}
|
|
|
|
pub(crate) fn sort_time<T: Copy + Ord>(xs: &[T], reverse: bool) -> Vec<usize> {
|
|
let mut x: Vec<(usize, T)> = xs.iter().enumerate().map(|(i, &t)| (i, t)).collect();
|
|
|
|
if reverse {
|
|
x.sort_by_key(|&(_, t)| Reverse(t));
|
|
} else {
|
|
x.sort_by_key(|&(_, t)| t);
|
|
}
|
|
|
|
x.into_iter().map(|(i, _)| i).collect()
|
|
}
|
|
|
|
/// Calculates the match quality of the given teams. A result is the draw probability in the association
|
|
///
|
|
/// Supports any number of groups. Values range roughly `[0, 1]`; 1 means a
|
|
/// perfectly balanced match.
|
|
///
|
|
/// # Panics
|
|
///
|
|
/// Panics if fewer than two teams are supplied, or if any group is
|
|
/// empty — match quality is a property of a contest between at least two
|
|
/// non-empty sides.
|
|
///
|
|
/// Also panics with "cannot invert a singular matrix" when every rating has
|
|
/// zero sigma *and* `beta` is zero. Nothing is then uncertain, so there is no
|
|
/// distribution to take the quality of; `Gaussian::from_ms(mu, 0.0)` is a point
|
|
/// mass and its `mu()` is not even well defined. Documented rather than
|
|
/// converted, because the input has no meaningful answer rather than an
|
|
/// awkward one.
|
|
#[must_use]
|
|
pub fn quality(teams: &[&[Gaussian]], beta: f64) -> f64 {
|
|
assert!(
|
|
teams.len() >= 2,
|
|
"quality() requires at least 2 teams, got {}",
|
|
teams.len()
|
|
);
|
|
assert!(
|
|
teams.iter().all(|group| !group.is_empty()),
|
|
"quality() requires every team to be non-empty"
|
|
);
|
|
|
|
let flatten_ratings = teams
|
|
.iter()
|
|
.flat_map(|group| group.iter())
|
|
.collect::<Vec<_>>();
|
|
|
|
let flatten_weights = vec![1.0; flatten_ratings.len()].into_boxed_slice();
|
|
|
|
let length = flatten_ratings.len();
|
|
|
|
let mut mean_matrix = Matrix::new(length, 1);
|
|
|
|
for (i, rating) in flatten_ratings.iter().enumerate() {
|
|
mean_matrix[(i, 0)] = rating.mu();
|
|
}
|
|
|
|
let mut variance_matrix = Matrix::new(length, length);
|
|
|
|
for (i, rating) in flatten_ratings.iter().enumerate() {
|
|
variance_matrix[(i, i)] = rating.sigma().powi(2);
|
|
}
|
|
|
|
let mut rotated_a_matrix = Matrix::new(teams.len() - 1, length);
|
|
|
|
// Row `row` contrasts group `row` (+weight) against group `row + 1`
|
|
// (-weight). `t` is the column where the current group's players start;
|
|
// the negative block begins immediately after it.
|
|
let mut t = 0;
|
|
|
|
for (row, group) in teams.windows(2).enumerate() {
|
|
let current = group[0];
|
|
let next = group[1];
|
|
|
|
for n in t..t + current.len() {
|
|
rotated_a_matrix[(row, n)] = flatten_weights[n];
|
|
}
|
|
|
|
t += current.len();
|
|
|
|
for n in t..t + next.len() {
|
|
rotated_a_matrix[(row, n)] = -flatten_weights[n];
|
|
}
|
|
}
|
|
|
|
let a_matrix = rotated_a_matrix.transpose();
|
|
|
|
let ata = beta.powi(2) * &rotated_a_matrix * &a_matrix;
|
|
let atsa = &rotated_a_matrix * &variance_matrix * &a_matrix;
|
|
|
|
let start = mean_matrix.transpose() * &a_matrix;
|
|
let middle = &ata + &atsa;
|
|
let end = &rotated_a_matrix * &mean_matrix;
|
|
|
|
let e_arg = (-0.5 * &start * &middle.inverse() * &end).determinant();
|
|
|
|
// `sqrt(det(ata) / det(middle))`, taken in log space. Both determinants are
|
|
// products of `k - 1` diagonal entries, so they leave `f64`'s range long
|
|
// before their ratio does: measured at the crate defaults, 150 groups was
|
|
// correct at `8.45e-53`, 200 returned `0`, and 250 returned `NaN` where the
|
|
// true value is `9.51e-88`. With a small beta it is sharper still — at
|
|
// `sigma = beta = 1e-3`, 60 groups returned `NaN` against a true `1.32e-9`.
|
|
//
|
|
// The ratio is what the answer needs and it is representable throughout, so
|
|
// the intermediates are the only thing that ever overflowed.
|
|
let ln_s_arg = ata.ln_abs_determinant() - middle.ln_abs_determinant();
|
|
|
|
libm::exp(e_arg + 0.5 * ln_s_arg)
|
|
}
|
|
|
|
#[cfg(test)]
|
|
mod tests {
|
|
/// The truncated variance must stay a variance across every branch, and
|
|
/// the branches must agree where they meet.
|
|
///
|
|
/// `v_w` now has three regimes for a tie — half-line, narrow-window, and
|
|
/// the direct subtraction — and a misplaced crossover between them is the
|
|
/// failure mode this guards. A jump at a boundary is visible here even
|
|
/// though the absolute values are not pinned.
|
|
#[test]
|
|
fn truncated_variance_is_continuous_across_the_tie_branches() {
|
|
for &alpha in &[50.0, 99.0, 100.0, 101.0, 1e3, 1e5, 1e6] {
|
|
// Sweep the window width across NARROW_WINDOW_RATIO and the
|
|
// half-line threshold, which sit at different widths per alpha.
|
|
let mut previous: Option<(f64, f64)> = None;
|
|
let mut width = alpha / (NARROW_WINDOW_RATIO * 100.0);
|
|
while width < 40.0 / alpha {
|
|
// mu = 0 puts the window at [-margin, margin]; shift it out to
|
|
// `alpha` by moving the mean instead.
|
|
let margin = width * 0.5;
|
|
let mu = -(alpha + width * 0.5);
|
|
let (v, one_minus_w) = v_w(mu, 1.0, margin, true);
|
|
|
|
assert!(v.is_finite(), "alpha {alpha}, width {width:e}: v = {v}");
|
|
assert!(
|
|
one_minus_w.is_finite() && one_minus_w > 0.0 && one_minus_w <= 1.0,
|
|
"alpha {alpha}, width {width:e}: 1 - w = {one_minus_w:e} is not a variance"
|
|
);
|
|
|
|
if let Some((pv, pw)) = previous {
|
|
// Consecutive widths differ by 2x, so the moments may not
|
|
// differ by more than a small multiple of that.
|
|
assert!(
|
|
one_minus_w / pw < 32.0 && pw / one_minus_w < 32.0,
|
|
"alpha {alpha}: 1 - w jumped from {pw:e} to {one_minus_w:e} \
|
|
at width {width:e} — a branch boundary is misplaced"
|
|
);
|
|
assert!(
|
|
(v - pv).abs() <= 8.0 * width.max(1e-12) + 1e-9 * v.abs(),
|
|
"alpha {alpha}: v jumped from {pv} to {v} at width {width:e}"
|
|
);
|
|
}
|
|
previous = Some((v, one_minus_w));
|
|
width *= 2.0;
|
|
}
|
|
}
|
|
}
|
|
|
|
/// The narrow-window form against high-precision quadrature.
|
|
///
|
|
/// These are the inputs where the direct `v^2 - u` subtraction had four
|
|
/// significant digits left and returned a negative variance.
|
|
#[test]
|
|
fn narrow_window_truncation_matches_quadrature() {
|
|
for &(alpha, d, expect_v, expect_w) in &[
|
|
(1e6, 2e-6, 1_000_000.000_000_687, 2.759_383_390_335_666e-13),
|
|
(1e4, 1e-6, 10_000.000_000_499_167, 8.333_291_666_831_727e-14),
|
|
(
|
|
1e3,
|
|
1e-5,
|
|
1_000.000_004_991_666_6,
|
|
8.333_291_666_803_818e-12,
|
|
),
|
|
] {
|
|
let (v, one_minus_w) = narrow_window_truncation(alpha, d);
|
|
assert!(
|
|
((v - expect_v) / expect_v).abs() < 1e-12,
|
|
"alpha {alpha:e}: v = {v}, want {expect_v}"
|
|
);
|
|
assert!(
|
|
((one_minus_w - expect_w) / expect_w).abs() < 1e-8,
|
|
"alpha {alpha:e}: 1 - w = {one_minus_w:e}, want {expect_w:e}"
|
|
);
|
|
}
|
|
}
|
|
|
|
/// A NaN must survive the fold from ANY position, not only the last.
|
|
///
|
|
/// The fold runs over a `HashMap`, so "last" is per-process hash order. The
|
|
/// end-to-end symptom was a NaN fit reported as `converged: true` in 16 of
|
|
/// 30 runs on identical input; these three cases are the deterministic form
|
|
/// of that, so a regression cannot hide behind a lucky seed.
|
|
#[test]
|
|
fn tuple_max_propagates_a_nan_from_any_position() {
|
|
let nan = (f64::NAN, f64::NAN);
|
|
let small = (1e-9, 1e-9);
|
|
let big = (1e-3, 1e-3);
|
|
|
|
// NaN last.
|
|
let step = tuple_max(tuple_max(big, small), nan);
|
|
assert!(!step_is_finite(step), "NaN last: {step:?}");
|
|
|
|
// NaN middle.
|
|
let step = tuple_max(tuple_max(big, nan), small);
|
|
assert!(!step_is_finite(step), "NaN middle: {step:?}");
|
|
|
|
// NaN first — the case a plain `>` comparison drops.
|
|
let step = tuple_max(tuple_max(nan, big), small);
|
|
assert!(!step_is_finite(step), "NaN first: {step:?}");
|
|
}
|
|
|
|
/// `f64::max` would pass the test above's first two cases and fail the
|
|
/// third, so pin that it is not what we use.
|
|
#[test]
|
|
fn tuple_max_is_not_f64_max() {
|
|
assert!(
|
|
f64::max(f64::NAN, 1.0) == 1.0,
|
|
"premise: f64::max drops NaN"
|
|
);
|
|
let (a, _) = tuple_max((f64::NAN, 0.0), (1.0, 0.0));
|
|
assert!(a.is_nan(), "tuple_max must not drop what f64::max drops");
|
|
}
|
|
|
|
/// Ordinary values are unaffected.
|
|
#[test]
|
|
fn tuple_max_still_takes_the_larger_component() {
|
|
assert_eq!(tuple_max((1.0, 5.0), (3.0, 2.0)), (3.0, 5.0));
|
|
assert_eq!(tuple_max((3.0, 2.0), (1.0, 5.0)), (3.0, 5.0));
|
|
}
|
|
|
|
use ::approx::assert_ulps_eq;
|
|
|
|
use super::*;
|
|
|
|
#[test]
|
|
fn test_sort_time() {
|
|
assert_eq!(sort_time(&[0i64, 1, 2, 0], true), vec![2, 1, 0, 3]);
|
|
}
|
|
|
|
/// Upper-tail values of the standard normal, from published tables. The
|
|
/// point is not the digits — these are 7-digit table values — but that a
|
|
/// number comes back at all: `1 - cdf` returned exactly zero for every one
|
|
/// of these.
|
|
#[test]
|
|
fn survival_function_survives_the_far_tail() {
|
|
for (z, expected) in [
|
|
(9.0f64, 1.128_588e-19),
|
|
(12.0, 1.776_482e-33),
|
|
(20.0, 2.753_624e-89),
|
|
(37.0, 5.725_571e-300),
|
|
] {
|
|
let got = sf(z, 0.0, 1.0);
|
|
assert!(got > 0.0, "sf({z}) collapsed to zero");
|
|
assert!(
|
|
(got - expected).abs() / expected < 1e-6, // published table values, 7 digits
|
|
"sf({z}) = {got}, expected ~{expected}"
|
|
);
|
|
assert_eq!(
|
|
1.0 - cdf(z, 0.0, 1.0),
|
|
0.0,
|
|
"the naive form should still be zero here"
|
|
);
|
|
}
|
|
}
|
|
|
|
/// Where no cancellation happens the two forms must agree exactly enough
|
|
/// that nothing else in the crate shifts.
|
|
#[test]
|
|
fn survival_function_matches_the_naive_form_where_that_form_works() {
|
|
for z in [-4.0f64, -1.0, 0.0, 0.5, 1.0, 2.0, 3.0, 4.0] {
|
|
let naive = 1.0 - cdf(z, 0.0, 1.0);
|
|
let direct = sf(z, 0.0, 1.0);
|
|
assert!(
|
|
(naive - direct).abs() < 1e-15,
|
|
"z={z}: naive {naive} vs direct {direct}"
|
|
);
|
|
}
|
|
}
|
|
|
|
#[test]
|
|
fn survival_and_cdf_partition_the_mass() {
|
|
for z in [-3.0f64, -0.5, 0.0, 1.0, 2.5] {
|
|
let total = sf(z, 1.0, 2.0) + cdf(z, 1.0, 2.0);
|
|
assert!((total - 1.0).abs() < 1e-15, "z={z}: {total}");
|
|
}
|
|
}
|
|
|
|
/// `erfcx` switches formulation at x = 2; the two sides must meet.
|
|
#[test]
|
|
fn erfcx_is_continuous_across_its_crossover() {
|
|
for x in [1.90f64, 1.99, 1.999, 2.0, 2.001, 2.01, 2.10] {
|
|
let direct = (x * x).exp() * erfc(x);
|
|
let scaled = erfcx(x);
|
|
assert!(
|
|
(direct - scaled).abs() / scaled < 1e-14,
|
|
"x={x}: direct {direct} vs erfcx {scaled}"
|
|
);
|
|
}
|
|
}
|
|
|
|
/// The whole reason `erfcx` exists: it stays finite and O(1/x) exactly
|
|
/// where `exp(x^2)` overflows and `erfc(x)` underflows.
|
|
#[test]
|
|
fn erfcx_stays_finite_where_its_factors_do_not() {
|
|
for x in [27.0f64, 50.0, 1.0e3, 1.0e8] {
|
|
let scaled = erfcx(x);
|
|
assert!(scaled.is_finite() && scaled > 0.0, "erfcx({x}) = {scaled}");
|
|
// Asymptotically erfcx(x) -> 1 / (x * sqrt(pi)).
|
|
let asymptote = 1.0 / (x * std::f64::consts::PI.sqrt());
|
|
assert!(
|
|
(scaled - asymptote).abs() / asymptote < 1e-2,
|
|
"erfcx({x}) = {scaled} strays from its asymptote {asymptote}"
|
|
);
|
|
assert!(
|
|
(x * x).exp().is_infinite(),
|
|
"x={x} should overflow the direct form"
|
|
);
|
|
}
|
|
}
|
|
|
|
/// Truncation must never produce a non-finite posterior. Before the scaled
|
|
/// formulation these returned NaN from `0 / 0` past about 39 sigma.
|
|
#[test]
|
|
fn truncation_stays_finite_arbitrarily_far_into_the_tail() {
|
|
for alpha in [0.0f64, 8.0, 38.0, 40.0, 100.0, 1.0e3, 1.0e6, 1.0e9, 1.0e15] {
|
|
for tie in [false, true] {
|
|
let (v, w) = v_w(-alpha, 1.0, if tie { 1.0 } else { 0.0 }, tie);
|
|
assert!(v.is_finite(), "alpha={alpha} tie={tie}: v = {v}");
|
|
assert!(w.is_finite(), "alpha={alpha} tie={tie}: w = {w}");
|
|
// sigma_trunc = sigma * sqrt(1 - w) must stay real.
|
|
assert!(
|
|
(0.0..=1.0).contains(&w),
|
|
"alpha={alpha} tie={tie}: w = {w} leaves sqrt(1 - w) imaginary"
|
|
);
|
|
|
|
let (mu_t, sigma_t) = trunc(-alpha, 1.0, if tie { 1.0 } else { 0.0 }, tie);
|
|
assert!(
|
|
mu_t.is_finite() && sigma_t.is_finite(),
|
|
"alpha={alpha} tie={tie}: trunc = ({mu_t}, {sigma_t})"
|
|
);
|
|
}
|
|
}
|
|
}
|
|
|
|
/// The Mills gap switches from subtraction to series at alpha = 100. Both
|
|
/// are supposed to be right there; if they disagree, the crossover is in
|
|
/// the wrong place.
|
|
#[test]
|
|
fn the_mills_gap_series_meets_the_scaled_form() {
|
|
for alpha in [50.0f64, 99.0, 100.0, 101.0, 200.0] {
|
|
let scaled = SQRT_2_OVER_PI / erfcx(alpha / SQRT_2) - alpha;
|
|
let inv = alpha.recip();
|
|
let inv_sq = inv * inv;
|
|
let series = inv * (1.0 - inv_sq * (2.0 - inv_sq * (10.0 - 74.0 * inv_sq)));
|
|
assert!(
|
|
(scaled - series).abs() / series < 1e-9,
|
|
"alpha={alpha}: scaled {scaled} vs series {series}"
|
|
);
|
|
}
|
|
}
|
|
|
|
/// Folding the tie branch to `mu <= 0` is only valid if v is odd in mu and
|
|
/// w is even. Assert the symmetry the implementation relies on.
|
|
#[test]
|
|
fn tie_truncation_is_odd_in_v_and_even_in_w() {
|
|
for mu in [0.5f64, 3.0, 20.0, 40.0, 100.0, 1.0e3] {
|
|
let (v_pos, w_pos) = v_w(mu, 1.0, 1.0, true);
|
|
let (v_neg, w_neg) = v_w(-mu, 1.0, 1.0, true);
|
|
assert!(
|
|
(v_pos + v_neg).abs() < 1e-9,
|
|
"mu={mu}: v should be odd, got {v_pos} and {v_neg}"
|
|
);
|
|
assert!(
|
|
(w_pos - w_neg).abs() < 1e-9,
|
|
"mu={mu}: w should be even, got {w_pos} and {w_neg}"
|
|
);
|
|
}
|
|
}
|
|
|
|
/// `erfc_inv`'s initial guess had the wrong sign, putting Newton on the
|
|
/// mirror image of the root. Three fixed iterations could not cross back,
|
|
/// so the error grew as the argument shrank: at `p_draw = 0.99` the margin
|
|
/// came out 0.503 where the answer is 2.576.
|
|
#[test]
|
|
fn erfc_inv_matches_known_quantiles() {
|
|
// sqrt(2) * erfc_inv(1 - p) is the standard normal quantile
|
|
// Phi^-1((1 + p) / 2).
|
|
for (p, exact) in [
|
|
(0.5f64, 0.674_489_750_196_081_7f64),
|
|
(0.9, 1.644_853_626_951_472_7),
|
|
(0.95, 1.959_963_984_540_054_2),
|
|
(0.99, 2.575_829_303_548_9),
|
|
(0.999, 3.290_526_731_491_896_4),
|
|
] {
|
|
let got = SQRT_2 * erfc_inv(1.0 - p);
|
|
assert!(
|
|
(got - exact).abs() / exact < 1e-14,
|
|
"p={p}: got {got}, exact {exact}"
|
|
);
|
|
}
|
|
}
|
|
|
|
/// The draw margin must grow with the draw probability. It did not: it ran
|
|
/// 0.674 -> 1.476 -> 0.503 -> 0.982 as `p_draw` went 0.5 -> 0.9 -> 0.99 ->
|
|
/// 0.999, which is not a rounding error but a broken function.
|
|
/// Deep in the tail the accuracy limit is the *caller's* argument, not this
|
|
/// function.
|
|
///
|
|
/// `compute_margin(0.999999, ..)` computes `1.0 - p_draw`, and 0.999999 is
|
|
/// not representable: the subtraction cancels and leaves 2.9e-11 of
|
|
/// relative error in the argument before `erfc_inv` is even entered. Given
|
|
/// an exactly-representable argument the result is good to 1.8e-16, so this
|
|
/// is inherent to taking `p_draw` near one rather than something to fix
|
|
/// here. At `p_draw = 0.999` the whole path is still accurate to 4e-16.
|
|
///
|
|
/// Worth pinning: measured against a 70-digit reference, `puruspe`'s
|
|
/// `inverfc` returns the identical wrong value for the identical reason,
|
|
/// which is what makes it clear the fault is upstream of both.
|
|
#[test]
|
|
fn erfc_inv_is_exact_given_an_exactly_representable_argument() {
|
|
// erfc(z / sqrt2) = 1e-6 exactly, so z = Phi^-1(0.9999995).
|
|
let got = SQRT_2 * erfc_inv(1e-6);
|
|
let exact = 4.891_638_475_698_59;
|
|
assert!(
|
|
(got - exact).abs() / exact < 1e-14,
|
|
"got {got}, exact {exact}"
|
|
);
|
|
}
|
|
|
|
#[test]
|
|
fn compute_margin_is_monotone_in_the_draw_probability() {
|
|
let mut previous = 0.0;
|
|
for p_draw in [
|
|
0.001f64, 0.01, 0.1, 0.25, 0.5, 0.75, 0.9, 0.99, 0.999, 0.9999,
|
|
] {
|
|
let margin = compute_margin(p_draw, 1.0);
|
|
assert!(
|
|
margin > previous,
|
|
"p_draw={p_draw}: margin {margin} did not exceed {previous}"
|
|
);
|
|
previous = margin;
|
|
}
|
|
}
|
|
|
|
/// Round-tripping the margin back through the model's own CDF must recover
|
|
/// the draw probability it was built from.
|
|
#[test]
|
|
fn compute_margin_round_trips_through_the_cdf() {
|
|
for p_draw in [0.001f64, 0.1, 0.5, 0.9, 0.99, 0.999] {
|
|
for sd in [0.5f64, 1.0, 5.892_557] {
|
|
let margin = compute_margin(p_draw, sd);
|
|
// P(|X| < margin) for X ~ N(0, sd^2).
|
|
let recovered = 1.0 - 2.0 * cdf(-margin, 0.0, sd);
|
|
assert!(
|
|
(recovered - p_draw).abs() < 1e-14,
|
|
"p_draw={p_draw} sd={sd}: recovered {recovered}"
|
|
);
|
|
}
|
|
}
|
|
}
|
|
|
|
/// `ln_pdf`, `ln_sf` and `ln_interval` exist so evidence stays exact where
|
|
/// the linear forms underflow. Past ~38 sigma the linear value is zero and
|
|
/// its log is whatever floor it was clamped to.
|
|
#[test]
|
|
fn log_space_helpers_stay_exact_where_the_linear_forms_underflow() {
|
|
for z in [40.0f64, 60.0, 100.0, 1000.0] {
|
|
assert_eq!(pdf(z, 0.0, 1.0), 0.0, "pdf should underflow at {z}");
|
|
assert_eq!(sf(z, 0.0, 1.0), 0.0, "sf should underflow at {z}");
|
|
|
|
let lp = ln_pdf(z, 0.0, 1.0);
|
|
let expected_lp = -(SQRT_TAU).ln() - 0.5 * z * z;
|
|
assert!(
|
|
(lp - expected_lp).abs() < 1e-9,
|
|
"ln_pdf({z}) = {lp}, expected {expected_lp}"
|
|
);
|
|
|
|
let ls = ln_sf(z, 0.0, 1.0);
|
|
// ln Phi(-z) ~ -z^2/2 - ln(z) - ln(sqrt(2 pi)) for large z.
|
|
let approx = -0.5 * z * z - z.ln() - SQRT_TAU.ln();
|
|
assert!(
|
|
(ls - approx).abs() / approx.abs() < 1e-3,
|
|
"ln_sf({z}) = {ls}, asymptote {approx}"
|
|
);
|
|
assert!(
|
|
ls < f64::MIN_POSITIVE.ln(),
|
|
"ln_sf({z}) still on the clamp floor"
|
|
);
|
|
}
|
|
}
|
|
|
|
/// Where nothing underflows, the log helpers must agree with the direct
|
|
/// forms exactly enough that nothing else in the crate shifts.
|
|
#[test]
|
|
fn log_space_helpers_agree_with_the_linear_forms_in_range() {
|
|
for z in [-3.0f64, -1.0, 0.0, 1.0, 2.0, 5.0, 10.0, 20.0] {
|
|
let lp = ln_pdf(z, 0.5, 2.0);
|
|
let direct_pdf = pdf(z, 0.5, 2.0);
|
|
assert!(
|
|
(lp.exp() - direct_pdf).abs() <= 1e-12 * direct_pdf,
|
|
"ln_pdf at {z}: {} vs {direct_pdf}",
|
|
lp.exp()
|
|
);
|
|
|
|
let ls = ln_sf(z, 0.5, 2.0);
|
|
let direct = sf(z, 0.5, 2.0);
|
|
assert!(
|
|
(ls.exp() - direct).abs() <= 1e-13 * direct.max(1e-300),
|
|
"ln_sf at {z}: {} vs {direct}",
|
|
ls.exp()
|
|
);
|
|
}
|
|
}
|
|
|
|
#[test]
|
|
fn ln_interval_matches_the_direct_difference_when_nothing_is_small() {
|
|
for mu in [-2.0f64, 0.0, 0.5, 2.0] {
|
|
let direct = cdf(1.0, mu, 1.0) - cdf(-1.0, mu, 1.0);
|
|
let logged = ln_interval(-1.0, 1.0, mu, 1.0).exp();
|
|
assert!(
|
|
(logged - direct).abs() <= 1e-13 * direct,
|
|
"mu={mu}: {logged} vs {direct}"
|
|
);
|
|
}
|
|
}
|
|
|
|
/// A window far out in the tail: both endpoints underflow together, so the
|
|
/// difference has to be taken in scaled form.
|
|
#[test]
|
|
fn ln_interval_survives_a_window_deep_in_the_tail() {
|
|
for mu in [-50.0f64, -100.0, -1000.0] {
|
|
let logged = ln_interval(-1.0, 1.0, mu, 1.0);
|
|
assert!(logged.is_finite(), "mu={mu}: {logged}");
|
|
assert!(
|
|
logged < f64::MIN_POSITIVE.ln(),
|
|
"mu={mu}: {logged} is stuck on the clamp floor"
|
|
);
|
|
// Dominated by the near edge: ln P ~ ln Phi(-(|mu| - 1)).
|
|
let near = ln_sf(-1.0, mu, 1.0);
|
|
assert!(
|
|
(logged - near).abs() < 5.0,
|
|
"mu={mu}: {logged} strays from the near-edge tail {near}"
|
|
);
|
|
}
|
|
}
|
|
|
|
#[test]
|
|
fn test_quality() {
|
|
let a = Gaussian::from_ms(25.0, 3.0);
|
|
let b = Gaussian::from_ms(25.0, 3.0);
|
|
|
|
let q = quality(&[&[a], &[b]], 25.0 / 3.0 / 2.0);
|
|
|
|
assert_ulps_eq!(q, 0.8115343414514944, epsilon = 1e-6)
|
|
}
|
|
}
|