logaritmiskandClaude Opus 5
327324c411
Twenty-one breaking changes in one release, with a live consumer. The
changelog lists them; `MIGRATING.md` says what to do about them, leading
with the three that change what an existing, *compiling* call returns —
unknown keys, predictions from a broken fit, and `Gaussian`'s operators
— since those are the ones the compiler will not find for you.
Every "after" snippet was compiled, not written from memory, and doing
so caught three errors in my own guide:
- `log_evidence_for(&[&"alice"])` does not compile at `K = String`. The
right spelling is `&["alice"]`, which works at *both* key types —
checked, because a guide that is right for half its readers is worse
than no guide.
- the same for `filtered_log_evidence_for`
- the `Analysis<'h> { joint: Joint<'h> }` example needs a history at the
default key type; pairing it with a `History<String>` does not compile
git-cliff skips merge commits now. Every branch lands with `--no-ff`, so
a release's merges outnumber its real commits and say nothing the merged
ones do not — 0.9.0's changelog had fourteen lines of them under "Other
(unconventional)". `ci:` commits get a group instead of falling through
to that catch-all.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011hcFjNDmHXZF8URGLku5zZ
1368 lines
54 KiB
Rust
1368 lines
54 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).
|
|
//!
|
|
//! Upgrading? `MIGRATING.md` in the repository root covers every breaking
|
|
//! change, with the ones that alter what an existing call *returns* called out
|
|
//! first.
|
|
//!
|
|
//! # 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 mod rating_rule;
|
|
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::{CompetitorField, InferenceError, OutcomeKind, Parameter, Shape, 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;
|
|
pub use rating_rule::{FnRule, NoRule, RatingRule, StartingPoint};
|
|
/// 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. Indices are assigned in interning order and
|
|
/// are stable for the life of a history; they are not portable between
|
|
/// histories, since the same key interns to a different slot under a different
|
|
/// ingestion order.
|
|
///
|
|
/// Crate-internal. It was public, along with `History::intern` and
|
|
/// `History::lookup` that produced one — and **nothing public ever accepted
|
|
/// one**, so it was a handle with nowhere to go. It also shadowed
|
|
/// `std::ops::Index`, which `CompetitorStore` implements. See #73.
|
|
#[derive(Copy, Clone, Default, PartialEq, PartialOrd, Eq, Ord, Hash, Debug)]
|
|
pub(crate) struct Index(usize);
|
|
|
|
impl Index {
|
|
/// The underlying slot number.
|
|
pub(crate) 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)
|
|
}
|
|
}
|