logaritmiskandClaude Opus 5
076a7ded8c
`Gaussian` publicly implemented `Mul`, `Div`, `Add` and `Sub`. They were
the EP product, cavity and variance-space convolutions, and every one of
them lies to a reader who takes the operator at face value:
a = N(10, 2) b = N(4, 3) c = N(1, 1)
a * b N(8.15, 1.66) not 40
a - b sigma GREW, 2 -> sqrt(4 + 9)
a * N(1, 0) mu = NaN "multiply by one"
a / c pi = -0.75 mu() prints a confident 0
The last is this crate's signature defect on a public operator. `Div` is
the cavity and can legitimately leave a negative precision, which is not
a distribution — and `mu()`/`sigma()` guard `pi <= 0` and report `0.0`
and `inf`, so it comes back as a plausible number with no panic, no
`Debug` marker and nothing to test against.
The four impls are now `pub(crate)` inherent methods that say what they
do: `ep_product`, `cavity`, `convolve`, `convolve_diff`, plus `scale`
for the one operation that genuinely is arithmetic. Nothing in a user's
workflow needed operator syntax; inference did, and it still has it.
`pi()` and `tau()` follow. Storing natural parameters is a performance
decision — it makes message passing two adds — not a contract. The
public surface is now exactly: `from_ms`, `from_mv`, `mu`, `sigma`,
`variance`, `probability_below`, `probability_above`. `from_mv` and
`variance` are promoted from `pub(crate)`; they are the honest pair for
callers who already hold a variance and should not pay a round trip
through the square root.
Four integration tests asserted bit-identity on `(pi, tau)`. They assert
it on `(mu, variance)` instead — still `assert_eq!`, still exact, and
`1/pi` and `tau/pi` are deterministic, so bit-equal natural parameters
give bit-equal moments. `a_nan_sigma_passes_through_from_ms` drops its
`|| g.pi().is_nan()` half: `sigma()` substitutes for `pi <= 0` and
`pi == inf`, so NaN survives to it only from a NaN precision.
`benches/gaussian.rs` is deleted. It timed two f64 additions through the
public operators, and keeping those public solely to feed it is the same
thing #73 objected to when a benchmark was dictating five public types.
The paths it covered are exercised by `batch` and `history_converge`
through the real call chain.
Closes #71.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011hcFjNDmHXZF8URGLku5zZ
580 lines
23 KiB
Rust
580 lines
23 KiB
Rust
use crate::{MU, N_INF, SIGMA};
|
|
|
|
/// A Gaussian distribution stored in natural parameters.
|
|
///
|
|
/// `pi = 1 / sigma^2` (precision)
|
|
/// `tau = mu * pi` (precision-adjusted mean)
|
|
///
|
|
/// Multiplication and division in message passing become pure adds/subs of
|
|
/// the stored fields with no `sqrt` or reciprocal in the hot path. `mu()` and
|
|
/// `sigma()` are accessors computed on demand.
|
|
#[derive(Clone, Copy, PartialEq, Debug)]
|
|
#[must_use]
|
|
pub struct Gaussian {
|
|
pi: f64,
|
|
tau: f64,
|
|
}
|
|
|
|
impl Gaussian {
|
|
/// Construct from mean and standard deviation.
|
|
///
|
|
/// # Panics
|
|
///
|
|
/// Panics if `sigma` is negative. NaN is deliberately allowed through: a
|
|
/// broken fit produces one, and `converge` reports that as
|
|
/// `NonFiniteResult` rather than panicking mid-inference.
|
|
///
|
|
/// A negative sigma used to be accepted and returned results **bit
|
|
/// identical** to its absolute value, because sigma only ever enters as
|
|
/// `sigma * sigma`. The sign was not rejected and not honoured; it simply
|
|
/// vanished. That is the same defect `HistoryBuilder::sigma`,
|
|
/// `HistoryBuilder::beta` and `Member::with_drift_scale` already reject.
|
|
///
|
|
/// # Very small sigma
|
|
///
|
|
/// `pi = 1 / sigma^2` leaves `f64`'s range below about `1.5e-154`, and
|
|
/// `tau = mu * pi` overflows sooner still — at a threshold that depends on
|
|
/// `mu`, so there is a band where `pi` is finite and only `tau` is not.
|
|
/// Both land on the same point-mass representation the `sigma == 0.0`
|
|
/// branch produces, and a point mass with a non-zero mean has `mu() = NaN`,
|
|
/// because `tau / pi` is `inf / inf`.
|
|
///
|
|
/// This is not rejected, because `approx` legitimately produces a very
|
|
/// small truncated sigma and inference must not panic. It is worth knowing
|
|
/// that such a `Gaussian` is not equal to itself, so two identical
|
|
/// declarations of one can be reported as conflicting.
|
|
pub const fn from_ms(mu: f64, sigma: f64) -> Self {
|
|
// NaN is admitted on purpose. A broken fit legitimately produces a NaN
|
|
// sigma — `sqrt` of a negative truncated variance — and the design is
|
|
// to propagate that to `converge`'s `NonFiniteResult` guard, not to
|
|
// panic inside inference. Rejecting it here turned that reporting path
|
|
// into a crash, which two tests caught immediately.
|
|
assert!(
|
|
sigma >= 0.0 || sigma.is_nan(),
|
|
"sigma must not be negative; it is only ever squared, so a negative \
|
|
value would silently behave as its absolute value"
|
|
);
|
|
if sigma == f64::INFINITY {
|
|
Self { pi: 0.0, tau: 0.0 }
|
|
} else if sigma == 0.0 {
|
|
// Point mass at mu. tau = mu * pi = mu * inf.
|
|
// For mu == 0 this is 0; for mu != 0 it is inf * mu = inf (IEEE).
|
|
// Only N00 (mu=0, sigma=0) is used in practice.
|
|
Self {
|
|
pi: f64::INFINITY,
|
|
tau: if mu == 0.0 { 0.0 } else { f64::INFINITY },
|
|
}
|
|
} else {
|
|
let pi = 1.0 / (sigma * sigma);
|
|
Self { pi, tau: mu * pi }
|
|
}
|
|
}
|
|
|
|
/// Construct from mean and *variance*, skipping the square-root round trip.
|
|
///
|
|
/// `from_ms(mu, var.sqrt())` immediately squares the root away again to
|
|
/// recover `pi = 1/var`. Variance-combining operations work in variance
|
|
/// space throughout, so they go through here instead and never take a
|
|
/// root. Use it whenever you already hold a variance —
|
|
/// [`variance`](Gaussian::variance) is its inverse.
|
|
#[inline]
|
|
pub fn from_mv(mu: f64, var: f64) -> Self {
|
|
if var == f64::INFINITY {
|
|
Self { pi: 0.0, tau: 0.0 }
|
|
} else if var == 0.0 {
|
|
// Point mass at mu; see `from_ms` for the tau convention.
|
|
Self {
|
|
pi: f64::INFINITY,
|
|
tau: if mu == 0.0 { 0.0 } else { f64::INFINITY },
|
|
}
|
|
} else {
|
|
let pi = 1.0 / var;
|
|
Self { pi, tau: mu * pi }
|
|
}
|
|
}
|
|
|
|
/// Construct directly from natural parameters.
|
|
#[inline]
|
|
pub(crate) const fn from_natural(pi: f64, tau: f64) -> Self {
|
|
Self { pi, tau }
|
|
}
|
|
|
|
/// Precision, `1 / sigma^2` — one of the two natural parameters.
|
|
///
|
|
/// This is the representation the type actually stores, which is why the EP
|
|
/// product and cavity (`Mul` / `Div`) are plain adds and subtracts. Larger
|
|
/// means more certain; `0.0` is an improper, uninformative message and
|
|
/// `inf` is a point mass.
|
|
#[inline]
|
|
pub(crate) fn pi(&self) -> f64 {
|
|
self.pi
|
|
}
|
|
|
|
/// Precision-adjusted mean, `mu / sigma^2` — the other natural parameter.
|
|
///
|
|
/// Stored rather than derived, for the same reason as [`Gaussian::pi`].
|
|
/// Meaningful only alongside `pi`: on its own it is not a location.
|
|
#[inline]
|
|
pub(crate) fn tau(&self) -> f64 {
|
|
self.tau
|
|
}
|
|
|
|
/// Mean skill: the point estimate.
|
|
///
|
|
/// Derived from the natural parameters as `tau / pi`. An improper message
|
|
/// (`pi <= 0`) has no defined mean and reports `0.0` — see
|
|
/// [`Gaussian::sigma`], which reports `inf` for the same state, and read
|
|
/// the two together before treating a mean as informative.
|
|
#[inline]
|
|
#[must_use]
|
|
pub fn mu(&self) -> f64 {
|
|
// A non-positive precision is an improper (uninformative) Gaussian — its mean is
|
|
// undefined. Treat it like `pi == 0` and return 0. EP message cancellation can land
|
|
// `pi` on a tiny negative value (round-off of exactly zero); without this guard
|
|
// `tau / pi` would yield a spurious finite mean.
|
|
if self.pi <= 0.0 {
|
|
0.0
|
|
} else {
|
|
self.tau / self.pi
|
|
}
|
|
}
|
|
|
|
/// Variance, without the root-and-square of `sigma().powi(2)`.
|
|
///
|
|
/// Mirrors [`sigma`](Gaussian::sigma)'s treatment of the improper
|
|
/// (infinite) and point-mass (zero) cases, and is the inverse of
|
|
/// [`from_mv`](Gaussian::from_mv).
|
|
#[inline]
|
|
#[must_use]
|
|
pub fn variance(&self) -> f64 {
|
|
if self.pi <= 0.0 {
|
|
f64::INFINITY
|
|
} else if self.pi.is_infinite() {
|
|
0.0
|
|
} else {
|
|
1.0 / self.pi
|
|
}
|
|
}
|
|
|
|
/// Standard deviation: how unsure this estimate is.
|
|
///
|
|
/// Derived as `1 / sqrt(pi)`. An improper message (`pi <= 0`) reports
|
|
/// `inf`, and a point mass (`pi == inf`) reports `0.0` — both are real
|
|
/// states rather than error codes, and both are legitimate for a converged
|
|
/// fit with degenerate parameters.
|
|
#[inline]
|
|
#[must_use]
|
|
pub fn sigma(&self) -> f64 {
|
|
// A non-positive precision is improper → infinite standard deviation. Guarding
|
|
// `pi <= 0.0` (not just `== 0.0`) keeps `1.0 / pi.sqrt()` from returning NaN when EP
|
|
// cancellation produces a tiny negative precision (round-off of exactly zero).
|
|
if self.pi <= 0.0 {
|
|
f64::INFINITY
|
|
} else if self.pi.is_infinite() {
|
|
0.0
|
|
} else {
|
|
1.0 / self.pi.sqrt()
|
|
}
|
|
}
|
|
|
|
/// How far this Gaussian moved from `other`, as `(|d mu|, |d sigma|)`.
|
|
///
|
|
/// Identical messages have not moved, whatever their parameters, and that
|
|
/// case is answered in natural space before touching `mu()`/`sigma()`. An
|
|
/// improper message has `pi == 0`, so `sigma()` is infinite — and
|
|
/// `inf - inf` is NaN, a NaN *change* for a message that did not change at
|
|
/// all. (`mu()` is guarded and returns 0.0 here, so the mean component was
|
|
/// never the problem; the sigma component alone produced `(0.0, NaN)`.)
|
|
///
|
|
/// That is reachable in ordinary inference: once a pairing is more than
|
|
/// about nine cavity-sigma apart the truncation is a no-op, `trunc / cavity`
|
|
/// is exactly the identity message, and the chain compares one identity
|
|
/// against another. Before this guard that produced `(0.0, NaN)`, which
|
|
/// silently disabled the sigma half of the convergence test.
|
|
pub(crate) fn delta(&self, other: Gaussian) -> (f64, f64) {
|
|
if self.pi == other.pi && self.tau == other.tau {
|
|
return (0.0, 0.0);
|
|
}
|
|
|
|
(
|
|
(self.mu() - other.mu()).abs(),
|
|
(self.sigma() - other.sigma()).abs(),
|
|
)
|
|
}
|
|
|
|
pub(crate) fn exclude(&self, other: Gaussian) -> Self {
|
|
let var = self.variance() - other.variance();
|
|
if var <= 0.0 {
|
|
// When sigma_self ≈ sigma_other (including ULP-level rounding differences
|
|
// from the pi→sigma accessor round-trip), the excluded contribution is N00.
|
|
// Computing from_ms(tiny_mu, 0.0) would give {pi:inf, tau:inf}, whose
|
|
// mu() = inf/inf = NaN. Returning N00 is correct: when both Gaussians
|
|
// carry the same variance, the residual is a point mass at 0.
|
|
return Gaussian::from_mv(0.0, 0.0);
|
|
}
|
|
|
|
Self::from_mv(self.mu() - other.mu(), var)
|
|
}
|
|
|
|
pub(crate) fn forget(&self, variance_delta: f64) -> Self {
|
|
Self::from_mv(self.mu(), self.variance() + variance_delta)
|
|
}
|
|
|
|
/// `P(X < x)` under this Gaussian.
|
|
///
|
|
/// The question a stopping rule asks: *how sure am I that this competitor's
|
|
/// true skill is below the cutoff?* Expressing that as a probability keeps
|
|
/// its meaning as sigma changes, where a `mu + z * sigma` band silently
|
|
/// means different confidence at different uncertainties — which is exactly
|
|
/// the regime a stopping rule operates in.
|
|
///
|
|
/// Accurate in the *lower* tail. For the upper tail use
|
|
/// [`Gaussian::probability_above`] rather than `1.0 - probability_below(x)`,
|
|
/// which cancels away every significant digit once the result is small.
|
|
///
|
|
/// An improper Gaussian (non-positive precision) has no defined mean, so
|
|
/// this returns `0.5` — the same convention `mu()` and `sigma()` follow.
|
|
#[must_use]
|
|
pub fn probability_below(&self, x: f64) -> f64 {
|
|
if self.pi <= 0.0 {
|
|
return 0.5;
|
|
}
|
|
crate::cdf(x, self.mu(), self.sigma())
|
|
}
|
|
|
|
/// `P(X > x)` under this Gaussian.
|
|
///
|
|
/// Computed as a survival function rather than `1 - cdf`, so it keeps full
|
|
/// relative precision in the upper tail: `1 - cdf` returns exactly zero
|
|
/// past about 8.3 sigma, where the true value is still 1e-19 and perfectly
|
|
/// representable. A stopping rule is evaluated precisely there — the
|
|
/// interesting cases are the ones near certainty.
|
|
///
|
|
/// An improper Gaussian returns `0.5`, as [`Gaussian::probability_below`].
|
|
#[must_use]
|
|
pub fn probability_above(&self, x: f64) -> f64 {
|
|
if self.pi <= 0.0 {
|
|
return 0.5;
|
|
}
|
|
crate::sf(x, self.mu(), self.sigma())
|
|
}
|
|
|
|
/// EP damping in natural-parameter space: `α·new + (1−α)·self`.
|
|
///
|
|
/// Used by within-game inference to stabilise oscillating fixed-point
|
|
/// loops on hard graphs. `alpha = 1.0` returns `new` exactly;
|
|
/// `alpha < 1.0` shrinks each per-step update.
|
|
pub(crate) fn damp_natural(self, new: Gaussian, alpha: f64) -> Gaussian {
|
|
Gaussian::from_natural(
|
|
alpha * new.pi() + (1.0 - alpha) * self.pi(),
|
|
alpha * new.tau() + (1.0 - alpha) * self.tau(),
|
|
)
|
|
}
|
|
}
|
|
|
|
impl Default for Gaussian {
|
|
fn default() -> Self {
|
|
Self::from_ms(MU, SIGMA)
|
|
}
|
|
}
|
|
|
|
impl Gaussian {
|
|
/// The EP factor **product**: multiply two messages about the same
|
|
/// variable.
|
|
///
|
|
/// Two natural-parameter additions and no square root, which is why the
|
|
/// type stores `pi` and `tau` rather than `mu` and `sigma`. This is the
|
|
/// hot path.
|
|
///
|
|
/// Not arithmetic — `N(10, 2).ep_product(N(4, 3))` is `N(8.15, 1.66)`,
|
|
/// nowhere near 40. It used to be spelled `a * b`, on a public `Mul` impl,
|
|
/// where that was a trap rather than a shorthand.
|
|
#[inline]
|
|
pub(crate) fn ep_product(self, rhs: Gaussian) -> Gaussian {
|
|
Self::from_natural(self.pi + rhs.pi, self.tau + rhs.tau)
|
|
}
|
|
|
|
/// The EP **cavity**: divide out a message this belief already absorbed.
|
|
///
|
|
/// The inverse of [`ep_product`](Gaussian::ep_product), and two
|
|
/// subtractions rather than two additions.
|
|
///
|
|
/// **May return an improper result.** Cancelling a message that carried
|
|
/// most of the precision leaves `pi <= 0`, which is not a distribution.
|
|
/// `mu()` reports `0.0` and `sigma()` reports `inf` for such a value —
|
|
/// both are the accessors' policy for "undefined", not answers. Measured:
|
|
/// `N(10, 2).cavity(N(1, 1))` has `pi = -0.75`, and its `mu()` prints a
|
|
/// confident `0`. That is why this is not a public operator.
|
|
#[inline]
|
|
pub(crate) fn cavity(self, rhs: Gaussian) -> Gaussian {
|
|
Self::from_natural(self.pi - rhs.pi, self.tau - rhs.tau)
|
|
}
|
|
|
|
/// Convolve two independent Gaussians: `N(mu1 + mu2, sqrt(v1 + v2))`.
|
|
///
|
|
/// The distribution of a *sum* of independent variables, so the variances
|
|
/// add — the result is always wider than either input. Used to combine a
|
|
/// skill with performance noise. Goes through `from_mv` and takes no root.
|
|
#[inline]
|
|
pub(crate) fn convolve(self, rhs: Gaussian) -> Gaussian {
|
|
Self::from_mv(self.mu() + rhs.mu(), self.variance() + rhs.variance())
|
|
}
|
|
|
|
/// Convolve a *difference*: `N(mu1 - mu2, sqrt(v1 + v2))`.
|
|
///
|
|
/// The means subtract and the variances still **add**, because a
|
|
/// difference of independent variables is no more certain than a sum. That
|
|
/// is the half that made the old `Sub` impl misleading: `a - b` grew the
|
|
/// sigma from 2 to `sqrt(4 + 9)`.
|
|
#[inline]
|
|
pub(crate) fn convolve_diff(self, rhs: Gaussian) -> Gaussian {
|
|
Self::from_mv(self.mu() - rhs.mu(), self.variance() + rhs.variance())
|
|
}
|
|
|
|
/// Scale by a constant: `mu` by `scalar`, `sigma` by `|scalar|`.
|
|
///
|
|
/// The one operation that *is* ordinary arithmetic — it is the
|
|
/// distribution of `scalar * X`. Used for per-member weights.
|
|
#[inline]
|
|
pub(crate) fn scale(self, scalar: f64) -> Gaussian {
|
|
if !scalar.is_finite() {
|
|
return N_INF;
|
|
}
|
|
if scalar == 0.0 {
|
|
// Scaling by 0 collapses to a point mass at 0 (sigma' = 0, mu' = 0).
|
|
// This is N00, the additive identity, NOT N_INF.
|
|
return Gaussian::from_mv(0.0, 0.0);
|
|
}
|
|
// sigma' = sigma * |scalar| => pi' = pi / scalar²
|
|
// mu' = mu * scalar => tau' = tau / scalar
|
|
Self::from_natural(self.pi / (scalar * scalar), self.tau / scalar)
|
|
}
|
|
}
|
|
|
|
#[cfg(test)]
|
|
mod tests {
|
|
/// A message that did not change must report no change, even when it is
|
|
/// improper. `mu()` of an improper Gaussian is `0/0 = NaN` and `sigma()` is
|
|
/// infinite, so the mean/sigma form reported `(NaN, NaN)` for two identical
|
|
/// identity messages — which silently disabled the sigma half of the
|
|
/// convergence test in `run_chain`.
|
|
#[test]
|
|
fn delta_of_two_identical_improper_messages_is_zero() {
|
|
let improper = crate::N_INF;
|
|
// `mu()` is guarded and returns 0.0 for an improper Gaussian, so the
|
|
// mean component was always fine. The NaN came from the sigma
|
|
// component alone: `inf - inf`. The pre-fix value was `(0.0, NaN)`.
|
|
assert!(improper.sigma().is_infinite(), "premise: sigma is infinite");
|
|
assert_eq!(improper.mu(), 0.0, "premise: mu is guarded, not NaN");
|
|
assert!(
|
|
(improper.sigma() - improper.sigma()).is_nan(),
|
|
"premise: the unguarded sigma difference is NaN"
|
|
);
|
|
assert_eq!(improper.delta(improper), (0.0, 0.0));
|
|
}
|
|
|
|
#[test]
|
|
fn delta_of_identical_proper_messages_is_zero() {
|
|
let g = Gaussian::from_ms(25.0, 8.0);
|
|
assert_eq!(g.delta(g), (0.0, 0.0));
|
|
}
|
|
|
|
/// The shortcut must not swallow a real difference.
|
|
#[test]
|
|
fn delta_still_measures_a_real_move() {
|
|
let a = Gaussian::from_ms(25.0, 8.0);
|
|
let b = Gaussian::from_ms(26.0, 9.0);
|
|
let (dmu, dsigma) = a.delta(b);
|
|
assert!((dmu - 1.0).abs() < 1e-12, "{dmu}");
|
|
assert!((dsigma - 1.0).abs() < 1e-12, "{dsigma}");
|
|
}
|
|
|
|
use super::*;
|
|
|
|
#[test]
|
|
fn non_positive_precision_is_improper_not_nan() {
|
|
// EP message cancellation can leave `pi` a tiny negative (round-off of exactly zero).
|
|
// Such a Gaussian is improper/uninformative: mu() must be 0 and sigma() infinite, not
|
|
// NaN. A NaN here propagates through the moment-space `Sub` in the game chain and
|
|
// poisons every skill in the slice.
|
|
let tiny_neg = Gaussian::from_natural(-5.55e-17, -8.88e-16);
|
|
assert_eq!(tiny_neg.mu(), 0.0);
|
|
assert!(tiny_neg.sigma().is_infinite());
|
|
|
|
// A frankly-negative precision is treated the same way.
|
|
let neg = Gaussian::from_natural(-1.0, 2.0);
|
|
assert_eq!(neg.mu(), 0.0);
|
|
assert!(neg.sigma().is_infinite());
|
|
|
|
// Subtracting such a message must not produce NaN (the original failure path).
|
|
let proper = Gaussian::from_ms(9.75, 1.256);
|
|
let diff = proper.convolve_diff(tiny_neg);
|
|
assert!(diff.pi().is_finite() && !diff.pi().is_nan());
|
|
assert!(diff.tau().is_finite() && !diff.tau().is_nan());
|
|
}
|
|
|
|
#[test]
|
|
fn convolve_adds_variances() {
|
|
let n = Gaussian::from_ms(25.0, 25.0 / 3.0);
|
|
let m = Gaussian::from_ms(0.0, 1.0);
|
|
let r = n.convolve(m);
|
|
assert!((r.mu() - 25.0).abs() < 1e-12);
|
|
assert!((r.sigma() - 8.393118874676116).abs() < 1e-10);
|
|
}
|
|
|
|
#[test]
|
|
fn convolve_diff_subtracts_means_and_adds_variances() {
|
|
let n = Gaussian::from_ms(25.0, 25.0 / 3.0);
|
|
let m = Gaussian::from_ms(1.0, 1.0);
|
|
let r = n.convolve_diff(m);
|
|
assert!((r.mu() - 24.0).abs() < 1e-12);
|
|
assert!((r.sigma() - 8.393118874676116).abs() < 1e-10);
|
|
}
|
|
|
|
#[test]
|
|
fn ep_product_is_not_arithmetic() {
|
|
let n = Gaussian::from_ms(25.0, 25.0 / 3.0);
|
|
let m = Gaussian::from_ms(0.0, 1.0);
|
|
let r = n.ep_product(m);
|
|
assert!((r.mu() - 0.35488958990536273).abs() < 1e-10);
|
|
assert!((r.sigma() - 0.992876838486922).abs() < 1e-10);
|
|
}
|
|
|
|
#[test]
|
|
fn cavity_undoes_a_product() {
|
|
let n = Gaussian::from_ms(25.0, 25.0 / 3.0);
|
|
let m = Gaussian::from_ms(0.0, 1.0);
|
|
let r = m.cavity(n);
|
|
assert!((r.mu() - (-0.3652597402597402)).abs() < 1e-10);
|
|
assert!((r.sigma() - 1.0072787050317253).abs() < 1e-10);
|
|
}
|
|
|
|
#[test]
|
|
fn test_n00_is_add_identity() {
|
|
// N00 (sigma=0) is the identity for `convolve`.
|
|
// N_INF (sigma=inf) is the identity for `ep_product`.
|
|
let g = Gaussian::from_ms(3.0, 2.0);
|
|
let n00 = Gaussian::from_ms(0.0, 0.0);
|
|
let r = n00.convolve(g);
|
|
assert!((r.mu() - g.mu()).abs() < 1e-12);
|
|
assert!((r.sigma() - g.sigma()).abs() < 1e-12);
|
|
}
|
|
|
|
#[test]
|
|
fn ep_product_adds_natural_parameters() {
|
|
// `ep_product` in nat-params should be pi_n + pi_m, tau_n + tau_m
|
|
let n = Gaussian::from_ms(2.0, 3.0);
|
|
let m = Gaussian::from_ms(1.0, 2.0);
|
|
let r = n.ep_product(m);
|
|
let expected_pi = n.pi() + m.pi();
|
|
let expected_tau = n.tau() + m.tau();
|
|
assert!((r.pi() - expected_pi).abs() < 1e-15);
|
|
assert!((r.tau() - expected_tau).abs() < 1e-15);
|
|
}
|
|
|
|
#[test]
|
|
fn cavity_subtracts_natural_parameters() {
|
|
let n = Gaussian::from_ms(2.0, 1.0);
|
|
let m = Gaussian::from_ms(1.0, 2.0);
|
|
let r = n.cavity(m);
|
|
let expected_pi = n.pi() - m.pi();
|
|
let expected_tau = n.tau() - m.tau();
|
|
assert!((r.pi() - expected_pi).abs() < 1e-15);
|
|
assert!((r.tau() - expected_tau).abs() < 1e-15);
|
|
}
|
|
|
|
#[test]
|
|
fn damp_natural_alpha_one_returns_new() {
|
|
let old = Gaussian::from_ms(1.0, 2.0);
|
|
let new = Gaussian::from_ms(5.0, 0.5);
|
|
let damped = old.damp_natural(new, 1.0);
|
|
assert_eq!(damped.pi(), new.pi());
|
|
assert_eq!(damped.tau(), new.tau());
|
|
}
|
|
|
|
#[test]
|
|
fn damp_natural_alpha_zero_returns_self() {
|
|
let old = Gaussian::from_ms(1.0, 2.0);
|
|
let new = Gaussian::from_ms(5.0, 0.5);
|
|
let damped = old.damp_natural(new, 0.0);
|
|
assert_eq!(damped.pi(), old.pi());
|
|
assert_eq!(damped.tau(), old.tau());
|
|
}
|
|
|
|
#[test]
|
|
fn damp_natural_alpha_half_is_midpoint_in_natural_params() {
|
|
let old = Gaussian::from_ms(1.0, 2.0);
|
|
let new = Gaussian::from_ms(5.0, 0.5);
|
|
let damped = old.damp_natural(new, 0.5);
|
|
let expected_pi = 0.5 * new.pi() + 0.5 * old.pi();
|
|
let expected_tau = 0.5 * new.tau() + 0.5 * old.tau();
|
|
assert!((damped.pi() - expected_pi).abs() < 1e-12);
|
|
assert!((damped.tau() - expected_tau).abs() < 1e-12);
|
|
}
|
|
}
|
|
|
|
#[cfg(test)]
|
|
mod tail_probability_tests {
|
|
use super::*;
|
|
|
|
#[test]
|
|
fn probability_below_matches_published_quantiles() {
|
|
let g = Gaussian::from_ms(0.0, 1.0);
|
|
for (x, expected) in [
|
|
(-1.959_963_984_540_054, 0.025),
|
|
(0.0, 0.5),
|
|
(1.281_551_565_544_6, 0.9),
|
|
(1.959_963_984_540_054, 0.975),
|
|
] {
|
|
let got = g.probability_below(x);
|
|
assert!(
|
|
(got - expected).abs() < 1e-12,
|
|
"P(X < {x}) = {got}, expected {expected}"
|
|
);
|
|
}
|
|
}
|
|
|
|
#[test]
|
|
fn the_two_tails_partition_the_mass() {
|
|
let g = Gaussian::from_ms(3.0, 2.0);
|
|
for x in [-4.0f64, 0.0, 3.0, 7.5] {
|
|
let total = g.probability_below(x) + g.probability_above(x);
|
|
assert!((total - 1.0).abs() < 1e-15, "at {x}: {total}");
|
|
}
|
|
}
|
|
|
|
/// The reason `probability_above` exists rather than `1 - probability_below`.
|
|
#[test]
|
|
fn probability_above_keeps_precision_where_the_complement_collapses() {
|
|
let g = Gaussian::from_ms(0.0, 1.0);
|
|
for (x, expected) in [(9.0f64, 1.128_588e-19), (20.0, 2.753_624e-89)] {
|
|
let got = g.probability_above(x);
|
|
assert!(
|
|
(got - expected).abs() / expected < 1e-6,
|
|
"P(X > {x}) = {got}, expected ~{expected}"
|
|
);
|
|
assert_eq!(
|
|
1.0 - g.probability_below(x),
|
|
0.0,
|
|
"the complement should still collapse at {x}"
|
|
);
|
|
}
|
|
}
|
|
|
|
#[test]
|
|
fn a_scaled_gaussian_shifts_and_stretches() {
|
|
let g = Gaussian::from_ms(25.0, 6.0);
|
|
assert!((g.probability_below(25.0) - 0.5).abs() < 1e-15);
|
|
// One sigma either side of the mean.
|
|
assert!((g.probability_below(31.0) - 0.841_344_746_068_543).abs() < 1e-12);
|
|
assert!((g.probability_above(19.0) - 0.841_344_746_068_543).abs() < 1e-12);
|
|
}
|
|
|
|
#[test]
|
|
fn an_improper_gaussian_is_uninformative_rather_than_nan() {
|
|
let improper = Gaussian::from_ms(0.0, f64::INFINITY);
|
|
assert_eq!(improper.probability_below(5.0), 0.5);
|
|
assert_eq!(improper.probability_above(5.0), 0.5);
|
|
}
|
|
}
|