logaritmiskandClaude Opus 5
8116fd081f
Scouting crates.io for a more accurate `erfc` than `libm` turned up a 1.16e-12 relative error in `erfc_inv` at `p_draw = 0.999999`, measured against a 70-digit `decimal` reference. It looked like too few Newton steps. It is not: adding a fourth changed nothing. The error is in forming the argument. `1.0 - 0.999999` is `1.0000000000287557e-06` — 0.999999 is not representable, and subtracting from one cancels, leaving 2.9e-11 of relative error before `erfc_inv` is entered. Given an exactly-representable argument it returns 1.8e-16. So the routine was never the problem, and the extra iteration has been reverted rather than shipped as a fix for a defect that was not there. `puruspe::inverfc` returns the identical wrong value for the identical reason, which is what makes the shared upstream cause obvious. Adds a test that separates the two, and corrects a quantile constant in `erfc_inv_matches_known_quantiles` that was recalled rather than computed: `Phi^-1(0.9999995)` is 4.89163847569859, not 4.891638475699099. The others were checked against the same reference and were right. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_011hcFjNDmHXZF8URGLku5zZ
1017 lines
37 KiB
Rust
1017 lines
37 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)]
|
|
|
|
/// 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},
|
|
};
|
|
|
|
#[cfg(feature = "approx")]
|
|
mod approx;
|
|
pub(crate) mod arena;
|
|
mod time;
|
|
mod time_slice;
|
|
pub use time_slice::{EventKind, TimeSlice};
|
|
mod acquisition;
|
|
mod color_group;
|
|
mod competitor;
|
|
mod convergence;
|
|
pub mod drift;
|
|
mod error;
|
|
mod event;
|
|
mod event_builder;
|
|
pub(crate) mod factor;
|
|
mod game;
|
|
pub mod gaussian;
|
|
pub mod graph;
|
|
mod history;
|
|
mod key_table;
|
|
mod matrix;
|
|
mod observer;
|
|
mod outcome;
|
|
mod predict;
|
|
pub(crate) mod quadrature;
|
|
mod rating;
|
|
pub(crate) mod schedule;
|
|
pub mod storage;
|
|
|
|
pub use acquisition::expected_information_gain;
|
|
pub use competitor::Competitor;
|
|
pub use convergence::{ConvergenceOptions, ConvergenceReport};
|
|
pub use drift::{ConstantDrift, Drift};
|
|
pub use error::InferenceError;
|
|
pub use event::{Event, Member, Team};
|
|
pub use event_builder::EventBuilder;
|
|
pub use game::{Game, GameOptions, OwnedGame};
|
|
pub use gaussian::Gaussian;
|
|
pub use history::{History, HistoryBuilder};
|
|
pub use key_table::KeyTable;
|
|
use matrix::Matrix;
|
|
pub use observer::{NullObserver, Observer};
|
|
pub use outcome::Outcome;
|
|
pub use predict::Prediction;
|
|
pub use rating::Rating;
|
|
pub use schedule::ScheduleReport;
|
|
pub use time::{Time, Untimed};
|
|
|
|
pub const BETA: f64 = 1.0;
|
|
pub const MU: f64 = 0.0;
|
|
pub const SIGMA: f64 = BETA * 6.0;
|
|
pub const GAMMA: f64 = BETA * 0.03;
|
|
pub const P_DRAW: f64 = 0.0;
|
|
pub const EPSILON: f64 = 1e-6;
|
|
pub const ITERATIONS: usize = 30;
|
|
|
|
/// Largest team count `History::predict_outcome` will enumerate.
|
|
///
|
|
/// The outcome space holds `n! * 2^(n-1)` events, so it grows factorially:
|
|
/// 1_920 at five teams, 23_040 at six, 322_560 at seven. Six is where
|
|
/// enumerating on a caller's behalf stops being reasonable.
|
|
pub const MAX_PREDICTED_TEAMS: usize = predict::MAX_TEAMS_FOR_DISTRIBUTION;
|
|
|
|
const SQRT_TAU: f64 = 2.5066282746310002;
|
|
/// `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.
|
|
const ASYMPTOTIC_MILLS_ALPHA: f64 = 100.0;
|
|
|
|
pub const N01: Gaussian = Gaussian::from_ms(0.0, 1.0);
|
|
pub const N00: Gaussian = Gaussian::from_ms(0.0, 0.0);
|
|
pub const N_INF: Gaussian = Gaussian::from_ms(0.0, f64::INFINITY);
|
|
|
|
#[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 gap = inv * (1.0 - inv_sq * (2.0 - inv_sq * (10.0 - 74.0 * inv_sq)));
|
|
let v = alpha + gap;
|
|
|
|
(v, v * gap)
|
|
}
|
|
|
|
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, 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.
|
|
if alpha >= ASYMPTOTIC_MILLS_ALPHA && alpha * (beta - alpha) >= HALF_LINE_WINDOW {
|
|
let (v, w) = half_line_truncation(alpha);
|
|
return (if flipped { -v } else { v }, 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,
|
|
)
|
|
};
|
|
|
|
let w = -(u - v.powi(2));
|
|
|
|
(if flipped { -v } else { v }, w)
|
|
}
|
|
}
|
|
|
|
fn trunc(mu: f64, sigma: f64, margin: f64, tie: bool) -> (f64, f64) {
|
|
let (v, w) = v_w(mu, sigma, margin, tie);
|
|
|
|
let mu_trunc = mu + sigma * v;
|
|
let sigma_trunc = sigma * (1.0 - 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)
|
|
}
|
|
|
|
pub(crate) fn tuple_max(v1: (f64, f64), v2: (f64, f64)) -> (f64, f64) {
|
|
(
|
|
if v1.0 > v2.0 { v1.0 } else { v2.0 },
|
|
if v1.1 > v2.1 { v1.1 } else { v2.1 },
|
|
)
|
|
}
|
|
|
|
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 rating groups. 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 rating groups are supplied, or if any group is
|
|
/// empty — match quality is a property of a contest between at least two
|
|
/// non-empty sides.
|
|
#[must_use]
|
|
pub fn quality(rating_groups: &[&[Gaussian]], beta: f64) -> f64 {
|
|
assert!(
|
|
rating_groups.len() >= 2,
|
|
"quality() requires at least 2 rating groups, got {}",
|
|
rating_groups.len()
|
|
);
|
|
assert!(
|
|
rating_groups.iter().all(|group| !group.is_empty()),
|
|
"quality() requires every rating group to be non-empty"
|
|
);
|
|
|
|
let flatten_ratings = rating_groups
|
|
.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(rating_groups.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 rating_groups.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();
|
|
let s_arg = ata.determinant() / middle.determinant();
|
|
|
|
libm::exp(e_arg) * s_arg.sqrt()
|
|
}
|
|
|
|
#[cfg(test)]
|
|
mod tests {
|
|
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)
|
|
}
|
|
}
|