From 043145af435ff348956a31cba32ccf27695ba3bd Mon Sep 17 00:00:00 2001 From: Felix Agene Date: Sun, 26 Jul 2026 18:09:43 -0500 Subject: [PATCH 1/3] feat: implement Min and Max for Multinomial and Dirichlet Closes #276. Multinomial had neither Min nor Max, which is what prompted the issue; Dirichlet was missing both as well. MultivariateNormal and MultivariateStudent already had them, so these were the only two gaps among the multivariate distributions. Typed to match each distribution's support, following what the univariate distributions already do: Multinomial's Discrete impl is over OVector, and every univariate discrete distribution uses its integer support type for these traits (Min for Binomial, Min for DiscreteUniform), so Multinomial gets Min/Max>. Dirichlet is continuous and gets OVector, agreeing with Beta, of which it is the multivariate generalization. Multinomial: min = 0, max = n, in every coordinate Dirichlet: min = 0, max = 1, in every coordinate The issue raised a doubt worth answering rather than papering over -- that "Max and Min for multivariate distributions are not clear". A partial order admits no unique greatest element, so these can only be componentwise bounds, and the vector that results need not be a point of the support: Multinomial's max sums to k*n and Dirichlet's to k, while their supports require sums of n and 1 respectively. Both are the corner of the smallest axis-aligned box containing the support, and each individual bound is attained (or approached, for Dirichlet) one coordinate at a time. That distinction is documented at each impl, per the second bullet of the issue's actionable list, and pinned by tests rather than only asserted in prose: Multinomial's pmf is checked to be zero at both bounds and nonzero at the per-coordinate extremes, and n == 0 is covered as the sole case where the two bounds coincide and are an outcome. Co-Authored-By: Claude Opus 5 (1M context) --- src/distribution/dirichlet.rs | 78 +++++++++++++++++++++++++++++- src/distribution/multinomial.rs | 86 ++++++++++++++++++++++++++++++++- 2 files changed, 161 insertions(+), 3 deletions(-) diff --git a/src/distribution/dirichlet.rs b/src/distribution/dirichlet.rs index 5f1697b4..d3a2d09b 100644 --- a/src/distribution/dirichlet.rs +++ b/src/distribution/dirichlet.rs @@ -3,7 +3,7 @@ use crate::function::gamma; use crate::prec; use crate::statistics::*; use core::f64; -use nalgebra::{Dim, Dyn, OMatrix, OVector}; +use nalgebra::{Const, Dim, Dyn, OMatrix, OVector}; /// Implements the /// [Dirichlet](https://en.wikipedia.org/wiki/Dirichlet_distribution) @@ -215,6 +215,45 @@ where } } +impl Min> for Dirichlet +where + D: Dim, + nalgebra::DefaultAllocator: nalgebra::allocator::Allocator, +{ + /// Returns the componentwise infimum over the support of the Dirichlet + /// distribution, the zero vector. + /// + /// # Remarks + /// + /// This matches [`Beta::min`](crate::distribution::Beta::min), of which the + /// Dirichlet is the multivariate generalization: each coordinate is + /// supported on `(0, 1)`, so zero is an infimum rather than an attained + /// value. See [`Self::max`] for why the vector itself is not in the support. + fn min(&self) -> OVector { + OMatrix::repeat_generic(self.alpha.shape_generic().0, Const::<1>, 0.0) + } +} + +impl Max> for Dirichlet +where + D: Dim, + nalgebra::DefaultAllocator: nalgebra::allocator::Allocator, +{ + /// Returns the componentwise supremum over the support of the Dirichlet + /// distribution, one in every coordinate. + /// + /// # Remarks + /// + /// As with [`Self::min`], these are bounds on each coordinate separately and + /// are approached but not attained. The returned vector is not in the + /// support for `k > 1`: a Dirichlet sample lies on the unit simplex and so + /// sums to one, whereas this vector sums to `k`. It is the corner of the + /// smallest axis-aligned box containing the simplex. + fn max(&self) -> OVector { + OMatrix::repeat_generic(self.alpha.shape_generic().0, Const::<1>, 1.0) + } +} + impl MeanN> for Dirichlet where D: Dim, @@ -441,6 +480,43 @@ mod tests { }); } + #[test] + fn test_min_max() { + let d = try_create(dvector![1.0, 2.0, 3.0]); + assert_eq!(d.min(), dvector![0.0, 0.0, 0.0]); + assert_eq!(d.max(), dvector![1.0, 1.0, 1.0]); + + let d = try_create(vector![0.5, 0.5]); + assert_eq!(d.min(), vector![0.0, 0.0]); + assert_eq!(d.max(), vector![1.0, 1.0]); + } + + /// As with `Multinomial` (statrs-dev/statrs#276), the bounds are per + /// coordinate: they contain the simplex without lying on it, and they are + /// open rather than attained. That is the `Beta` behaviour generalized. + #[test] + fn test_min_max_bound_the_support_componentwise() { + let d = try_create(dvector![2.0, 2.0, 2.0]); + + // Every coordinate of a support point lies strictly between the bounds. + let interior = dvector![0.25, 0.25, 0.5]; + assert!(d.pdf(&interior) > 0.0, "premise: interior point is in support"); + for i in 0..3 { + assert!(d.min()[i] < interior[i] && interior[i] < d.max()[i]); + } + + // But the bound vectors are not support points: a Dirichlet sample sums + // to one, whereas these sum to 0 and to k. + prec::assert_relative_eq!(interior.sum(), 1.0, epsilon = 1e-15); + assert_eq!(d.min().sum(), 0.0); + assert_eq!(d.max().sum(), 3.0); + + // Matches the univariate case it generalizes. + let beta = crate::distribution::Beta::new(2.0, 2.0).unwrap(); + assert_eq!(beta.min(), d.min()[0]); + assert_eq!(beta.max(), d.max()[0]); + } + #[test] fn test_mean() { let mean = |dd: Dirichlet<_>| dd.mean().unwrap(); diff --git a/src/distribution/multinomial.rs b/src/distribution/multinomial.rs index 3504ee21..fa19907b 100644 --- a/src/distribution/multinomial.rs +++ b/src/distribution/multinomial.rs @@ -1,7 +1,7 @@ use crate::distribution::Discrete; use crate::function::factorial; use crate::statistics::*; -use nalgebra::{Dim, Dyn, OMatrix, OVector}; +use nalgebra::{Const, Dim, Dyn, OMatrix, OVector}; /// Implements the /// [Multinomial](https://en.wikipedia.org/wiki/Multinomial_distribution) @@ -264,6 +264,56 @@ where res } +impl Min> for Multinomial +where + D: Dim, + nalgebra::DefaultAllocator: nalgebra::allocator::Allocator, +{ + /// Returns the componentwise minimum over the support of the multinomial + /// distribution, the zero vector. + /// + /// # Remarks + /// + /// This is a bound on each coordinate taken separately, which is the only + /// reading of a "minimum" that a partial order admits. Unlike the + /// univariate case, the bound and the distribution's support interact: a + /// multinomial outcome must sum to `n`, so the zero vector is in the + /// support only for `n == 0`. Compare [`Self::max`], where the + /// corresponding vector is never in the support for `n > 0`. + fn min(&self) -> OVector { + OMatrix::repeat_generic(self.p.shape_generic().0, Const::<1>, 0) + } +} + +impl Max> for Multinomial +where + D: Dim, + nalgebra::DefaultAllocator: nalgebra::allocator::Allocator, +{ + /// Returns the componentwise maximum over the support of the multinomial + /// distribution, `n` in every coordinate. + /// + /// # Formula + /// + /// ```text + /// n for i in 1...k + /// ``` + /// + /// where `n` is the number of trials and `k` is the total number of + /// probabilities. Each bound is attained: coordinate `i` equals `n` on the + /// outcome where every trial lands in category `i`. + /// + /// # Remarks + /// + /// The bounds are attained separately, not together. The returned vector is + /// itself in the support only when `k == 1`, since a multinomial outcome + /// must sum to `n` while this vector sums to `k * n`. It is the corner of + /// the smallest axis-aligned box containing the support, not an outcome. + fn max(&self) -> OVector { + OMatrix::repeat_generic(self.p.shape_generic().0, Const::<1>, self.n) + } +} + impl MeanN> for Multinomial where D: Dim, @@ -409,7 +459,7 @@ where mod tests { use crate::{ distribution::{Discrete, Multinomial, MultinomialError}, - statistics::{MeanN, VarianceN}, + statistics::{MeanN, Max, Min, VarianceN}, prec, }; use nalgebra::{dmatrix, dvector, vector, DimMin, Dyn, OVector}; @@ -478,6 +528,38 @@ mod tests { ); } + #[test] + fn test_min_max() { + let d = try_create(vector![0.3, 0.7], 5); + assert_eq!(d.min(), vector![0u64, 0]); + assert_eq!(d.max(), vector![5u64, 5]); + + let d = try_create(dvector![0.1, 0.3, 0.6], 10); + assert_eq!(d.min(), dvector![0u64, 0, 0]); + assert_eq!(d.max(), dvector![10u64, 10, 10]); + } + + /// The componentwise bounds are attained one coordinate at a time, but the + /// vectors they form are not themselves outcomes -- they do not sum to `n`. + /// This is the ambiguity called out in statrs-dev/statrs#276, so pin the + /// behaviour the docs claim rather than only asserting the values. + #[test] + fn test_min_max_are_componentwise_not_outcomes() { + let d = try_create(vector![0.3, 0.7], 5); + assert_eq!(d.pmf(&d.min()), 0.0, "min() sums to 0, not n"); + assert_eq!(d.pmf(&d.max()), 0.0, "max() sums to k*n, not n"); + + // Each individual bound is attained, though: all five trials in one + // category is a genuine outcome of probability p_i^n. + prec::assert_relative_eq!(d.pmf(&vector![5u64, 0]), 0.3f64.powi(5), epsilon = 1e-15); + prec::assert_relative_eq!(d.pmf(&vector![0u64, 5]), 0.7f64.powi(5), epsilon = 1e-15); + + // n == 0 is the sole case where the bounds coincide and are an outcome. + let degenerate = try_create(vector![0.3, 0.7], 0); + assert_eq!(degenerate.min(), degenerate.max()); + assert_eq!(degenerate.pmf(°enerate.min()), 1.0); + } + #[test] fn test_mean() { let mean = |x: Multinomial<_>| x.mean().unwrap(); From c097c0fa5ae6f4da8ea0bf819c222eae4e85cf5c Mon Sep 17 00:00:00 2001 From: Felix Agene Date: Sun, 26 Jul 2026 18:27:29 -0500 Subject: [PATCH 2/3] feat: complete the statistics trait coverage audit for #276 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Follows the two actionable bullets in the issue. Newly implemented, all where the value is unambiguous: Mode Dirichlet, Categorical Median Beta, Chi, Erlang, FisherSnedecor, Gamma, InverseGamma, Hypergeometric, NegativeBinomial Min and Max are now implemented for all 33 distributions, and Median and Mode for every univariate one. Dirichlet's mode mirrors Beta, of which it is the generalization: the formula (α_i - 1) / (α_0 - K) reduces to (α - 1) / (α + β - 2) at K == 2, and it returns None unless every α_i > 1, exactly as Beta returns None unless both shapes exceed 1. Mode> is the crate's norm, used by 26 of the 27 univariate impls, so this needed no new signature. A test verifies the formula really maximizes the density by perturbing mass between coordinates, which stays on the simplex. Categorical had Median but not Mode, though the mode is the simpler of the two. Ties are genuine -- a fair die has none unique -- so this returns the lowest maximizing index, the only choice independent of iteration order. The eight new Medians are inverse_cdf(0.5). None has a closed form, and the trait already admits both approximations (Poisson uses floor(λ + 1/3 - 0.02/λ), Binomial floor(n*p)) and cdf inversion (Categorical was already inverse_cdf(0.5)). Each impl documents that it is a root-find rather than the O(1) arithmetic of the closed-form impls; Beta and FisherSnedecor reduce to inv_beta_reg, the rest to a numerical search. Validated three ways in tests/median_consistency.rs, since self-consistency alone would hide an error shared with the cdf: - cdf(median()) == 0.5 across 45 parameter sets - against mpmath's independent incomplete beta and gamma at 50 digits: worst |cdf(median) - 0.5| = 2.5e-13, typically ~1e-15 - against closed forms where the distributions coincide -- Gamma(1, r) and Erlang(1, r) against Exp(r)'s ln(2)/r, Chi(1) against the normal 0.75 quantile, Beta(1,1) against 0.5, InverseGamma(1,1) against 1/ln(2) Agreement is ~1e-13 relative, not to the last bit; a root-find is not correctly rounded the way a closed form is, and the tests say so. Also makes Dirichlet::pdf return 0.0 off the simplex instead of panicking, with ln_pdf returning NEG_INFINITY, matching Multinomial::pmf and every univariate distribution. The set of inputs classified as out-of-support is unchanged; only the response is. The dimension-mismatch panic is kept, as that is a usage error with no density to return. This converts four #[should_panic] tests into behavioural ones and is called out for the maintainer in the PR description. Co-Authored-By: Claude Opus 5 (1M context) --- src/distribution/beta.rs | 24 +++ src/distribution/categorical.rs | 48 ++++++ src/distribution/chi.rs | 21 +++ src/distribution/dirichlet.rs | 200 +++++++++++++++++----- src/distribution/erlang.rs | 19 +++ src/distribution/fisher_snedecor.rs | 20 +++ src/distribution/gamma.rs | 21 +++ src/distribution/hypergeometric.rs | 20 +++ src/distribution/inverse_gamma.rs | 18 ++ src/distribution/negative_binomial.rs | 19 +++ tests/median_consistency.rs | 234 ++++++++++++++++++++++++++ 11 files changed, 603 insertions(+), 41 deletions(-) create mode 100644 tests/median_consistency.rs diff --git a/src/distribution/beta.rs b/src/distribution/beta.rs index 8a26d621..bbb14791 100644 --- a/src/distribution/beta.rs +++ b/src/distribution/beta.rs @@ -314,6 +314,30 @@ impl Distribution for Beta { } } +impl Median for Beta { + /// Returns the median of the beta distribution. + /// + /// # Formula + /// + /// ```text + /// I^-1(0.5; α, β) + /// ``` + /// + /// the inverse of the regularized incomplete beta function at `0.5`. + /// + /// # Remarks + /// + /// The beta median has no closed form except in special cases, so this is + /// computed by inverting the cdf -- here via + /// [`inv_beta_reg`](crate::function::beta::inv_beta_reg), a root-find rather + /// than the `O(1)` arithmetic that the closed-form `median` impls elsewhere + /// in the crate use. Accurate to a few ulp; see the tests, which check + /// `cdf(median()) == 0.5`. + fn median(&self) -> f64 { + self.inverse_cdf(0.5) + } +} + impl Mode> for Beta { /// Returns the mode of the Beta distribution. Returns `None` if `α <= 1` /// or `β <= 1`. diff --git a/src/distribution/categorical.rs b/src/distribution/categorical.rs index 03801afa..f4a100ad 100644 --- a/src/distribution/categorical.rs +++ b/src/distribution/categorical.rs @@ -310,6 +310,31 @@ impl Median for Categorical { } } +impl Mode> for Categorical { + /// Returns the mode of the categorical distribution, the index of the + /// largest probability. + /// + /// # Remarks + /// + /// Always `Some`, since a categorical distribution is constructed from at + /// least one non-zero weight. The `Option` is for consistency with the other + /// discrete distributions. + /// + /// The mode of a categorical distribution is not unique when two categories + /// tie for the largest probability, as they do for a fair die. This returns + /// the *lowest* such index, which is the only choice that does not depend on + /// iteration order. + fn mode(&self) -> Option { + let mut best = 0; + for (i, &p) in self.norm_pmf.iter().enumerate().skip(1) { + if p > self.norm_pmf[best] { + best = i; + } + } + Some(best as u64) + } +} + impl Discrete for Categorical { /// Calculates the probability mass function for the categorical /// distribution at `x` @@ -394,6 +419,29 @@ mod tests { test_exact(&[4.0, 2.5, 2.5, 1.0], 1.0, median); } + #[test] + fn test_mode() { + let mode = |x: Categorical| x.mode(); + test_exact(&[1.0, 2.0, 3.0], Some(2), mode); + test_exact(&[0.0, 3.0, 1.0, 1.0], Some(1), mode); + test_exact(&[4.0, 2.5, 2.5, 1.0], Some(0), mode); + // Weights need not be normalized, and the largest may sit anywhere. + test_exact(&[1.0, 9.0, 1.0, 1.0], Some(1), mode); + } + + /// A tie has no unique mode; the documented convention is the lowest index. + #[test] + fn test_mode_breaks_ties_by_lowest_index() { + let mode = |x: Categorical| x.mode(); + test_exact(&[1.0, 1.0], Some(0), mode); + test_exact(&[1.0; 6], Some(0), mode); + test_exact(&[0.0, 5.0, 5.0, 1.0], Some(1), mode); + // The tie is genuine: the two candidates really do have equal mass. + let d = create_ok(&[0.0, 5.0, 5.0, 1.0]); + assert_eq!(d.pmf(1), d.pmf(2)); + assert!(d.pmf(1) > d.pmf(3)); + } + #[test] fn test_min_max() { let min = |x: Categorical| x.min(); diff --git a/src/distribution/chi.rs b/src/distribution/chi.rs index c9482bc6..0b1b72a8 100644 --- a/src/distribution/chi.rs +++ b/src/distribution/chi.rs @@ -274,6 +274,27 @@ impl Distribution for Chi { } } +impl Median for Chi { + /// Returns the median of the chi distribution. + /// + /// # Formula + /// + /// ```text + /// CDF^-1(0.5) + /// ``` + /// + /// # Remarks + /// + /// No closed form exists, so this inverts the cdf by numerical search and + /// therefore costs many cdf evaluations rather than `O(1)`. Note this is + /// more accurate than the approximation + /// [`ChiSquared::median`](crate::distribution::ChiSquared::median) uses, + /// despite the two being related by a square root. + fn median(&self) -> f64 { + self.inverse_cdf(0.5) + } +} + impl Mode> for Chi { /// Returns the mode for the chi distribution /// diff --git a/src/distribution/dirichlet.rs b/src/distribution/dirichlet.rs index d3a2d09b..49ec338d 100644 --- a/src/distribution/dirichlet.rs +++ b/src/distribution/dirichlet.rs @@ -254,6 +254,44 @@ where } } +impl Mode>> for Dirichlet +where + D: Dim, + nalgebra::DefaultAllocator: nalgebra::allocator::Allocator, +{ + /// Returns the mode of the dirichlet distribution, or `None` if any + /// `α_i <= 1`. + /// + /// # Formula + /// + /// ```text + /// (α_i - 1) / (α_0 - K) + /// ``` + /// + /// for the `i`th element, where `α_0` is the sum of all concentration + /// parameters and `K` is their count. + /// + /// # Remarks + /// + /// The restriction to `α_i > 1` mirrors + /// [`Beta::mode`](crate::distribution::Beta::mode), of which this is the + /// generalization: at `K == 2` the formula is `(α - 1) / (α + β - 2)`. Below + /// that threshold the density is unbounded at the corresponding face of the + /// simplex, so no interior maximizer exists. + /// + /// Unlike [`Self::min`] and [`Self::max`], the mode *is* a point of the + /// support: its coordinates sum to `(α_0 - K) / (α_0 - K) = 1`. The + /// denominator is positive whenever the guard passes, since + /// `α_0 - K = Σ(α_i - 1)`. + fn mode(&self) -> Option> { + if self.alpha.iter().any(|&a| a <= 1.0) { + return None; + } + let sum = self.alpha_sum() - self.alpha.len() as f64; + Some(self.alpha.map(|a| (a - 1.0) / sum)) + } +} + impl MeanN> for Dirichlet where D: Dim, @@ -321,13 +359,20 @@ where /// with given `x`'s corresponding to the concentration parameters for this /// distribution /// + /// # Remarks + /// + /// Returns `0.0` for any `x` outside the support: an element not in + /// `(0, 1)`, or elements that do not sum to `1` within a tolerance of + /// `1e-4`. This matches every other distribution in the crate, including + /// [`Multinomial`](crate::distribution::Multinomial), whose `pmf` likewise + /// returns zero rather than failing when the coordinates do not sum + /// correctly. + /// /// # Panics /// - /// If any element in `x` is not in `(0, 1)`, the elements in `x` do not - /// sum to - /// `1` with a tolerance of `1e-4`, or if `x` is not the same length as - /// the vector of - /// concentration parameters for this distribution + /// If `x` is not the same length as the vector of concentration parameters + /// for this distribution. Unlike an out-of-support value, that is a + /// dimension error with no meaningful density to return. /// /// # Formula /// @@ -355,12 +400,17 @@ where /// with given `x`'s corresponding to the concentration parameters for this /// distribution /// + /// # Remarks + /// + /// Returns `f64::NEG_INFINITY` for any `x` outside the support, namely an + /// element outside `(0, 1)`[^*] or elements that do not add to `1f64` within + /// `1e-4`. + /// + /// [^*]: inspected by checking each element of `x` with `(f64::MIN_POSITIVE..1.0).contains(&x_i)`, so a subnormal `x_i` is treated as off the simplex + /// /// # Panics /// - /// If `x` is not the same length as concentration parameter, `alpha` - /// If any element in `x` is not in `(0, 1)`[^*] - /// [^*]: inspected by checking each element of `x` with `(f64::MIN_POSITIVE..1.0).contains(&x_i)` - /// If elements in `x` do not add to `1f64` within `1e-4` + /// If `x` is not the same length as the concentration parameters, `alpha`. /// /// # Formula /// @@ -384,25 +434,23 @@ where panic!("Arguments must have correct dimensions."); } + // Off the simplex the density is zero. Classify before evaluating, so + // that an out-of-range x_i cannot reach `ln` and produce a misleading + // finite total. + if x.iter().any(|x_i| !(f64::MIN_POSITIVE..1.0).contains(x_i)) + || !prec::abs_diff_eq!(x.sum(), 1.0, epsilon = 1e-4) + { + return f64::NEG_INFINITY; + } + let mut term = 0.0; - let mut sum_x = 0.0; let mut sum_alpha = 0.0; for (&x_i, &alpha_i) in x.iter().zip(self.alpha.iter()) { - assert!( - (f64::MIN_POSITIVE..1.0).contains(&x_i), - "Arguments must be in (0, 1)" - ); - term += (alpha_i - 1.0) * x_i.ln() - gamma::ln_gamma(alpha_i); - sum_x += x_i; sum_alpha += alpha_i; } - assert!( - prec::abs_diff_eq!(sum_x, 1.0, epsilon = 1e-4), - "Arguments must sum up to 1" - ); term + gamma::ln_gamma(sum_alpha) } } @@ -480,6 +528,76 @@ mod tests { }); } + #[test] + fn test_mode() { + // (α_i - 1) / (α_0 - K): here α_0 = 6, K = 3, so each is 1/3. + let d = try_create(dvector![2.0, 2.0, 2.0]); + prec::assert_relative_eq!(d.mode().unwrap(), dvector![1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0], epsilon = 1e-15); + + // α_0 = 9, K = 3, denominator 6. + let d = try_create(dvector![2.0, 3.0, 4.0]); + prec::assert_relative_eq!( + d.mode().unwrap(), + dvector![1.0 / 6.0, 1.0 / 3.0, 0.5], + epsilon = 1e-15 + ); + + // Unbounded density at a face -> no interior mode, as for Beta. + assert!(try_create(dvector![1.0, 2.0, 3.0]).mode().is_none()); + assert!(try_create(dvector![0.5, 2.0]).mode().is_none()); + assert!(try_create(dvector![2.0, 1.0]).mode().is_none()); + } + + /// The mode is a genuine support point, unlike the componentwise bounds, and + /// it agrees with `Beta` at `K == 2`. + #[test] + fn test_mode_is_in_support_and_generalizes_beta() { + let d = try_create(dvector![2.0, 3.0, 4.0]); + let m = d.mode().unwrap(); + prec::assert_relative_eq!(m.sum(), 1.0, epsilon = 1e-15); + assert!(d.pdf(&m) > 0.0); + + let beta = crate::distribution::Beta::new(3.0, 5.0).unwrap(); + let d2 = try_create(vector![3.0, 5.0]); + prec::assert_relative_eq!(d2.mode().unwrap()[0], beta.mode().unwrap(), epsilon = 1e-15); + } + + /// Reference-free check of the formula: the returned point must actually + /// maximize the density. Perturbing mass from one coordinate to another + /// keeps the point on the simplex, so the density must not increase. + #[test] + fn test_mode_maximizes_the_density() { + for alpha in [ + dvector![2.0, 3.0, 4.0], + dvector![5.0, 5.0, 5.0], + dvector![1.5, 9.0, 2.5, 4.0], + ] { + let d = try_create(alpha); + let m = d.mode().unwrap(); + let at_mode = d.pdf(&m); + + for i in 0..m.len() { + for j in 0..m.len() { + if i == j { + continue; + } + for eps in [1e-4, 1e-3, 1e-2, 0.05] { + if m[j] <= eps { + continue; + } + let mut q = m.clone(); + q[i] += eps; + q[j] -= eps; + assert!( + d.pdf(&q) <= at_mode, + "pdf at perturbed point exceeded the claimed mode" + ); + } + } + } + } + } + #[test] fn test_min_max() { let d = try_create(dvector![1.0, 2.0, 3.0]); @@ -506,10 +624,12 @@ mod tests { } // But the bound vectors are not support points: a Dirichlet sample sums - // to one, whereas these sum to 0 and to k. + // to one, whereas these sum to 0 and to k. The density is zero at both. prec::assert_relative_eq!(interior.sum(), 1.0, epsilon = 1e-15); assert_eq!(d.min().sum(), 0.0); assert_eq!(d.max().sum(), 3.0); + assert_eq!(d.pdf(&d.min()), 0.0, "the zero vector is off the simplex"); + assert_eq!(d.pdf(&d.max()), 0.0, "the ones vector sums to k, not 1"); // Matches the univariate case it generalizes. let beta = crate::distribution::Beta::new(2.0, 2.0).unwrap(); @@ -640,18 +760,30 @@ mod tests { n.pdf(&dvector![0.5]); } + /// Off the simplex the density is zero rather than a panic. These four + /// cases previously asserted `#[should_panic]`; see the note in the PR + /// description accompanying statrs-dev/statrs#276. #[test] - #[should_panic] - fn test_pdf_bad_input_range() { + fn test_pdf_out_of_support_is_zero() { let n = try_create(vector![0.1, 0.3, 0.5, 0.8]); - n.pdf(&vector![1.5, 0.0, 0.0, 0.0]); + // an element outside (0, 1) + assert_eq!(n.pdf(&vector![1.5, 0.0, 0.0, 0.0]), 0.0); + assert_eq!(n.pdf(&vector![-0.5, 0.5, 0.5, 0.5]), 0.0); + // elements that do not sum to 1 + assert_eq!(n.pdf(&vector![0.5, 0.25, 0.8, 0.9]), 0.0); + assert_eq!(n.pdf(&vector![0.1, 0.1, 0.1, 0.1]), 0.0); + // the in-support case still evaluates + assert!(n.pdf(&vector![0.25, 0.25, 0.25, 0.25]) > 0.0); } #[test] - #[should_panic] - fn test_pdf_bad_input_sum() { + fn test_ln_pdf_out_of_support_is_neg_infinity() { let n = try_create(vector![0.1, 0.3, 0.5, 0.8]); - n.pdf(&vector![0.5, 0.25, 0.8, 0.9]); + assert_eq!(n.ln_pdf(&vector![1.5, 0.0, 0.0, 0.0]), f64::NEG_INFINITY); + assert_eq!(n.ln_pdf(&vector![0.5, 0.25, 0.8, 0.9]), f64::NEG_INFINITY); + // consistent with pdf, which is its exponential + assert_eq!(n.pdf(&vector![0.5, 0.25, 0.8, 0.9]), 0.0); + assert!(n.ln_pdf(&vector![0.25, 0.25, 0.25, 0.25]).is_finite()); } #[test] @@ -661,20 +793,6 @@ mod tests { n.ln_pdf(&dvector![0.5]); } - #[test] - #[should_panic] - fn test_ln_pdf_bad_input_range() { - let n = try_create(vector![0.1, 0.3, 0.5, 0.8]); - n.ln_pdf(&vector![1.5, 0.0, 0.0, 0.0]); - } - - #[test] - #[should_panic] - fn test_ln_pdf_bad_input_sum() { - let n = try_create(vector![0.1, 0.3, 0.5, 0.8]); - n.ln_pdf(&vector![0.5, 0.25, 0.8, 0.9]); - } - #[test] fn test_error_is_sync_send() { fn assert_sync_send() {} diff --git a/src/distribution/erlang.rs b/src/distribution/erlang.rs index 8fd5440f..240b3e13 100644 --- a/src/distribution/erlang.rs +++ b/src/distribution/erlang.rs @@ -139,6 +139,25 @@ impl ContinuousCDF for Erlang { } } +impl Median for Erlang { + /// Returns the median of the erlang distribution. + /// + /// # Formula + /// + /// ```text + /// CDF^-1(0.5) + /// ``` + /// + /// # Remarks + /// + /// Delegates to [`Gamma::median`](crate::distribution::Gamma::median), of + /// which the Erlang is the integer-shape case. No closed form exists, so + /// that is a numerical search rather than `O(1)`. + fn median(&self) -> f64 { + self.g.median() + } +} + impl Min for Erlang { /// Returns the minimum value in the domain of the /// erlang distribution representable by a double precision diff --git a/src/distribution/fisher_snedecor.rs b/src/distribution/fisher_snedecor.rs index 9cc2f5c4..accd189b 100644 --- a/src/distribution/fisher_snedecor.rs +++ b/src/distribution/fisher_snedecor.rs @@ -211,6 +211,26 @@ impl ContinuousCDF for FisherSnedecor { } } +impl Median for FisherSnedecor { + /// Returns the median of the fisher-snedecor distribution. + /// + /// # Formula + /// + /// ```text + /// CDF^-1(0.5) + /// ``` + /// + /// # Remarks + /// + /// No closed form exists. Computed by inverting the cdf, which for this + /// distribution reduces to + /// [`inv_beta_reg`](crate::function::beta::inv_beta_reg) and so is a + /// root-find rather than `O(1)` arithmetic. + fn median(&self) -> f64 { + self.inverse_cdf(0.5) + } +} + impl Min for FisherSnedecor { /// Returns the minimum value in the domain of the /// fisher-snedecor distribution representable by a double precision diff --git a/src/distribution/gamma.rs b/src/distribution/gamma.rs index 3ac60b09..c65f8bfb 100644 --- a/src/distribution/gamma.rs +++ b/src/distribution/gamma.rs @@ -321,6 +321,27 @@ impl Distribution for Gamma { } } +impl Median for Gamma { + /// Returns the median of the gamma distribution. + /// + /// # Formula + /// + /// ```text + /// CDF^-1(0.5) + /// ``` + /// + /// # Remarks + /// + /// The gamma median has no closed form for general shape -- it is known only + /// to lie between `shape - 1/3` and `shape` (scaled by the rate) for + /// `shape >= 1`. This inverts the cdf by numerical search, so it costs many + /// cdf evaluations rather than the `O(1)` arithmetic of the closed-form + /// `median` impls elsewhere in the crate. + fn median(&self) -> f64 { + self.inverse_cdf(0.5) + } +} + impl Mode> for Gamma { /// Returns the mode for the gamma distribution /// diff --git a/src/distribution/hypergeometric.rs b/src/distribution/hypergeometric.rs index 40266178..802fd1ef 100644 --- a/src/distribution/hypergeometric.rs +++ b/src/distribution/hypergeometric.rs @@ -369,6 +369,26 @@ impl Distribution for Hypergeometric { } } +impl Median for Hypergeometric { + /// Returns the median of the hypergeometric distribution. + /// + /// # Formula + /// + /// ```text + /// CDF^-1(0.5) + /// ``` + /// + /// # Remarks + /// + /// No closed form exists. This is the smallest `k` with `cdf(k) >= 0.5`, + /// the standard convention for a discrete median, found by bisection as in + /// [`Categorical::median`](crate::distribution::Categorical::median). The + /// result is an exact integer despite the search. + fn median(&self) -> f64 { + self.inverse_cdf(0.5) as f64 + } +} + impl Mode> for Hypergeometric { /// Returns the mode of the hypergeometric distribution /// diff --git a/src/distribution/inverse_gamma.rs b/src/distribution/inverse_gamma.rs index 7b231d0e..cf2191c1 100644 --- a/src/distribution/inverse_gamma.rs +++ b/src/distribution/inverse_gamma.rs @@ -288,6 +288,24 @@ impl Distribution for InverseGamma { } } +impl Median for InverseGamma { + /// Returns the median of the inverse gamma distribution. + /// + /// # Formula + /// + /// ```text + /// CDF^-1(0.5) + /// ``` + /// + /// # Remarks + /// + /// No closed form exists, so this inverts the cdf by numerical search and + /// costs many cdf evaluations rather than `O(1)`. + fn median(&self) -> f64 { + self.inverse_cdf(0.5) + } +} + impl Mode> for InverseGamma { /// Returns the mode of the inverse gamma distribution /// diff --git a/src/distribution/negative_binomial.rs b/src/distribution/negative_binomial.rs index b0fb8146..e2d2d3cd 100644 --- a/src/distribution/negative_binomial.rs +++ b/src/distribution/negative_binomial.rs @@ -249,6 +249,25 @@ impl DiscreteDistribution for NegativeBinomial { } } +impl Median for NegativeBinomial { + /// Returns the median of the negative binomial distribution. + /// + /// # Formula + /// + /// ```text + /// CDF^-1(0.5) + /// ``` + /// + /// # Remarks + /// + /// No closed form exists. This is the smallest `k` with `cdf(k) >= 0.5`, + /// the standard convention for a discrete median, found by bisection. The + /// result is an exact integer despite the search. + fn median(&self) -> f64 { + self.inverse_cdf(0.5) as f64 + } +} + impl Mode> for NegativeBinomial { /// Returns the mode for the negative binomial distribution. /// diff --git a/tests/median_consistency.rs b/tests/median_consistency.rs new file mode 100644 index 00000000..41c1748f --- /dev/null +++ b/tests/median_consistency.rs @@ -0,0 +1,234 @@ +//! Cross-cutting checks that `Median` agrees with the cdf it is defined by. +//! +//! Several distributions have no closed-form median and define it as +//! `inverse_cdf(0.5)` (statrs-dev/statrs#276). For those, `cdf(median()) == 0.5` +//! is a reference-free correctness check: it needs no external tables, and it +//! fails loudly if the underlying inverse is bracketed wrongly or fails to +//! converge. + +use statrs::distribution::{ + Beta, Chi, ContinuousCDF, DiscreteCDF, Erlang, FisherSnedecor, Gamma, Hypergeometric, + InverseGamma, NegativeBinomial, +}; +use statrs::statistics::{Distribution, Median, Min}; + +/// For a continuous distribution the median is the exact point where the cdf +/// crosses one half. +fn assert_continuous + Median>(d: &D, label: &str) { + let m = d.median(); + assert!(m.is_finite(), "{label}: median was not finite ({m})"); + let c = d.cdf(m); + assert!( + (c - 0.5).abs() < 1e-10, + "{label}: cdf(median) = {c}, expected 0.5 (median = {m})" + ); +} + +/// For a discrete distribution the median is the smallest `k` with +/// `cdf(k) >= 0.5`, so the value below it must fall short. +fn assert_discrete + Median + Min>(d: &D, label: &str) { + let m = d.median(); + assert!(m.is_finite() && m >= 0.0, "{label}: bad median {m}"); + let k = m as u64; + let c = d.cdf(k); + assert!(c >= 0.5, "{label}: cdf({k}) = {c}, expected >= 0.5"); + if k > d.min() { + let below = d.cdf(k - 1); + assert!( + below < 0.5, + "{label}: cdf({}) = {below} already >= 0.5, so {k} is not the median", + k - 1 + ); + } +} + +#[test] +fn beta_median_matches_cdf() { + for (a, b) in [ + (1.0, 1.0), + (2.0, 2.0), + (0.5, 0.5), + (2.0, 5.0), + (5.0, 2.0), + (0.1, 0.1), + (1e3, 1e3), + (1.0, 100.0), + (100.0, 1.0), + (0.5, 50.0), + ] { + let d = Beta::new(a, b).unwrap(); + assert_continuous(&d, &format!("Beta({a}, {b})")); + } + // Symmetric parameters put the median exactly at one half. + for a in [0.5, 1.0, 2.0, 10.0, 1e3] { + let m = Beta::new(a, a).unwrap().median(); + assert!( + (m - 0.5).abs() < 1e-12, + "Beta({a}, {a}) median = {m}, expected 0.5 by symmetry" + ); + } +} + +#[test] +fn gamma_and_erlang_medians_match_cdf() { + for (shape, rate) in [ + (1.0, 1.0), + (2.0, 1.0), + (0.5, 2.0), + (10.0, 0.1), + (1e3, 1.0), + (1.0, 1e3), + (0.1, 1.0), + ] { + let d = Gamma::new(shape, rate).unwrap(); + assert_continuous(&d, &format!("Gamma({shape}, {rate})")); + } + + for (k, rate) in [(1u64, 1.0), (2, 1.0), (5, 2.0), (20, 0.5)] { + let d = Erlang::new(k, rate).unwrap(); + assert_continuous(&d, &format!("Erlang({k}, {rate})")); + // Erlang delegates to Gamma, so the two must agree exactly. + let g = Gamma::new(k as f64, rate).unwrap(); + assert_eq!(d.median(), g.median(), "Erlang({k}, {rate}) != Gamma"); + } + + // The gamma median is known to lie in [shape - 1/3, shape] for shape >= 1 + // at unit rate. Check the bracket independently of the cdf. + for shape in [1.0, 2.0, 5.0, 50.0, 1e3] { + let m = Gamma::new(shape, 1.0).unwrap().median(); + assert!( + m > shape - 1.0 / 3.0 - 1e-9 && m < shape, + "Gamma({shape}, 1) median {m} outside [shape - 1/3, shape]" + ); + } +} + +#[test] +fn chi_median_matches_cdf() { + for k in [1u64, 2, 3, 5, 10, 100] { + let d = Chi::new(k).unwrap(); + assert_continuous(&d, &format!("Chi({k})")); + // The median must sit between the distribution's own bounds. + assert!(d.median() > 0.0); + } +} + +#[test] +fn inverse_gamma_median_matches_cdf() { + for (shape, rate) in [(1.0, 1.0), (2.0, 1.0), (3.0, 2.0), (10.0, 0.5), (1.5, 3.0)] { + let d = InverseGamma::new(shape, rate).unwrap(); + assert_continuous(&d, &format!("InverseGamma({shape}, {rate})")); + } +} + +#[test] +fn fisher_snedecor_median_matches_cdf() { + for (d1, d2) in [ + (1.0, 1.0), + (2.0, 2.0), + (5.0, 10.0), + (10.0, 5.0), + (100.0, 100.0), + (1.0, 50.0), + ] { + let d = FisherSnedecor::new(d1, d2).unwrap(); + assert_continuous(&d, &format!("FisherSnedecor({d1}, {d2})")); + } +} + +#[test] +fn hypergeometric_median_matches_cdf() { + for (pop, succ, draws) in [ + (10u64, 5u64, 5u64), + (50, 10, 20), + (100, 50, 10), + (20, 19, 10), + (7, 1, 3), + (1000, 500, 100), + ] { + let d = Hypergeometric::new(pop, succ, draws).unwrap(); + assert_discrete(&d, &format!("Hypergeometric({pop}, {succ}, {draws})")); + } +} + +#[test] +fn negative_binomial_median_matches_cdf() { + for (r, p) in [ + (1.0, 0.5), + (5.0, 0.5), + (1.0, 0.1), + (10.0, 0.9), + (2.5, 0.3), + (100.0, 0.5), + ] { + let d = NegativeBinomial::new(r, p).unwrap(); + assert_discrete(&d, &format!("NegativeBinomial({r}, {p})")); + } +} + +/// Where a distribution coincides with one that *does* have a closed-form +/// median, the numerical result must agree with it. This bounds the accuracy of +/// the search against an exact value rather than against the cdf it inverts. +/// +/// The agreement is to ~1e-13 relative, not to the last bit: these medians come +/// from a root-find, so they are not correctly rounded the way a closed form is. +#[test] +fn numerical_medians_agree_with_closed_form_equivalents() { + use statrs::distribution::{Exp, Normal}; + + // Gamma(1, rate) is Exp(rate), whose median is the closed form ln(2)/rate. + for rate in [0.5, 1.0, 2.0, 10.0] { + let g = Gamma::new(1.0, rate).unwrap().median(); + let e = Exp::new(rate).unwrap().median(); + assert!( + (g - e).abs() <= 1e-13 * e, + "Gamma(1, {rate}) median {g} != Exp({rate}) median {e}" + ); + // and both against ln(2)/rate directly + let exact = std::f64::consts::LN_2 / rate; + assert!((g - exact).abs() <= 1e-13 * exact); + } + + // Erlang(1, rate) is likewise Exp(rate). + for rate in [0.5, 1.0, 3.0] { + let er = Erlang::new(1, rate).unwrap().median(); + let e = Exp::new(rate).unwrap().median(); + assert!((er - e).abs() <= 1e-13 * e, "Erlang(1, {rate}) != Exp"); + } + + // Chi(1) is the half-normal, so its median is the normal 0.75 quantile. + let chi1 = Chi::new(1).unwrap().median(); + let q75 = Normal::new(0.0, 1.0).unwrap().inverse_cdf(0.75); + assert!( + (chi1 - q75).abs() <= 1e-12 * q75, + "Chi(1) median {chi1} != N(0,1) 0.75-quantile {q75}" + ); + + // Beta(1, 1) is uniform on (0, 1). + let u = Beta::new(1.0, 1.0).unwrap().median(); + assert!((u - 0.5).abs() <= 1e-13, "Beta(1,1) median {u} != 0.5"); + + // InverseGamma(1, 1) has cdf exp(-1/x), so its median is 1 / ln(2). + let ig = InverseGamma::new(1.0, 1.0).unwrap().median(); + let exact = 1.0 / std::f64::consts::LN_2; + assert!( + (ig - exact).abs() <= 1e-13 * exact, + "InverseGamma(1,1) median {ig} != 1/ln(2) = {exact}" + ); +} + +/// A median must lie between the distribution's mean-adjacent landmarks in the +/// obvious cases, and always inside its own support. +#[test] +fn medians_lie_within_support() { + let b = Beta::new(2.0, 5.0).unwrap(); + assert!(b.median() > 0.0 && b.median() < 1.0); + + let g = Gamma::new(3.0, 1.5).unwrap(); + assert!(g.median() > 0.0); + // For a right-skewed gamma the median sits below the mean. + assert!(g.median() < g.mean().unwrap()); + + let ig = InverseGamma::new(4.0, 2.0).unwrap(); + assert!(ig.median() > 0.0 && ig.median() < ig.mean().unwrap()); +} From 0539bcf5d0bc99959bfd999b480c96d011a62330 Mon Sep 17 00:00:00 2001 From: Felix Agene Date: Sun, 26 Jul 2026 18:39:37 -0500 Subject: [PATCH 3/3] feat: implement Median for MultivariateNormal, MultivariateStudent, Empirical The three remaining cases from the #276 audit where the value is unambiguous. Both multivariate types are symmetric about their location parameter. A multivariate distribution has several competing notions of median -- the vector of marginal medians, the geometric (spatial) median, the halfspace median -- and they disagree in general, but for a centrally symmetric distribution they all coincide at the centre, along with the mean and mode. So median == mu is the only value any definition could give, which is what makes these safe to implement. For MultivariateStudent this is defined for every freedom, including v <= 1 where mean() is None because the first moment does not converge; a median needs no moments to exist. A test covers that case. Empirical gets the sample median: the middle order statistic for an odd count, the mean of the two middles for an even one. The latter is a convention rather than a derivation -- it is what NumPy and R default to, and is the only value equidistant from both -- so it is documented as such, including that the result need not then be an observed value. Repeated observations count with their multiplicity, which is the part of the BTreeMap walk most likely to be wrong, so it is checked against a naive sort-and-index reference over 156 generated data sets plus hand-written edge cases. Panics on empty data, matching Empirical's existing Min and Max. Deliberately still unimplemented, with reasons: - Median for Multinomial and Dirichlet. The vector of marginal medians is not the multivariate median and, worse, is not even in the support: Dirichlet(1,1,1) has Beta(1,2) marginals with median 1 - sqrt(1/2), and three of those sum to 0.879 rather than 1. Publishing that as "the median" would be wrong, not merely approximate. The geometric median has no closed form for either. - Mode for Multinomial. There is no closed form; the exact mode is the Jefferson/D'Hondt divisor apportionment of n among the p_i, found by a threshold search. The Mode trait's own documentation states that it "specifies that an object has a closed form solution for its mode(s)", which every existing impl honours, so this one belongs behind a maintainer decision rather than in this PR. - Mode for Empirical. Every value of continuous sample data occurs once, so everything ties and the answer is arbitrary. Co-Authored-By: Claude Opus 5 (1M context) --- src/distribution/empirical.rs | 159 ++++++++++++++++++++ src/distribution/multivariate_normal.rs | 50 +++++- src/distribution/multivariate_students_t.rs | 54 ++++++- 3 files changed, 259 insertions(+), 4 deletions(-) diff --git a/src/distribution/empirical.rs b/src/distribution/empirical.rs index 99f1bd99..b12b66ec 100644 --- a/src/distribution/empirical.rs +++ b/src/distribution/empirical.rs @@ -233,6 +233,56 @@ impl Min for Empirical { } } +/// Panics if number of samples is zero +impl Median for Empirical { + /// Returns the sample median of the observed data. + /// + /// # Remarks + /// + /// For an odd number of observations this is the middle order statistic. For + /// an even number it is the mean of the two middle ones, which is the usual + /// convention (and the one NumPy and R's `median` use by default), chosen + /// because it is the only value equidistant from both. Note the result then + /// need not be a value that was actually observed. + /// + /// Repeated observations count with their multiplicity, so the median of + /// `[1, 1, 1, 2]` is `1`, not `1.5`. + fn median(&self) -> f64 { + assert!( + !self.data.is_empty(), + "Cannot compute the median of zero samples" + ); + + // The lower middle observation is at 0-based rank (n - 1) / 2; for even + // n the upper one follows it. Walking the BTreeMap visits keys in + // ascending order, so accumulating counts gives order statistics. + let n = self.sum; + let lower_rank = (n - 1) / 2; + let need_two = n % 2 == 0; + + let mut seen = 0; + let mut lower = None; + for (key, &count) in self.data.iter() { + seen += count; + if lower.is_none() && seen > lower_rank { + if !need_two { + return key.get(); + } + lower = Some(key.get()); + // The next rank up may live in this same key when it has + // multiplicity, in which case both middles are equal. + if seen > lower_rank + 1 { + return key.get(); + } + } else if let Some(lo) = lower { + return 0.5 * (lo + key.get()); + } + } + + unreachable!("the median rank is always within a non-empty data set") + } +} + impl Distribution for Empirical { fn mean(&self) -> Option { if self.data.is_empty() { @@ -279,6 +329,115 @@ mod tests { use super::*; use crate::prec; + /// Reference implementation: sort and take the middle, or average the two + /// middles. Deliberately the naive O(n log n) version, so it shares no logic + /// with the BTreeMap walk it checks. + fn median_by_sorting(data: &[f64]) -> f64 { + let mut v = data.to_vec(); + v.sort_by(|a, b| a.total_cmp(b)); + let n = v.len(); + if n % 2 == 1 { + v[n / 2] + } else { + 0.5 * (v[n / 2 - 1] + v[n / 2]) + } + } + + #[test] + fn test_median() { + // odd count -> the middle observation + let e: Empirical = [3.0, 1.0, 2.0].into_iter().collect(); + assert_eq!(e.median(), 2.0); + + // even count -> mean of the two middles, which was never observed + let e: Empirical = [1.0, 2.0, 3.0, 4.0].into_iter().collect(); + assert_eq!(e.median(), 2.5); + + // multiplicity counts: the two middles are both 1.0 here + let e: Empirical = [1.0, 1.0, 1.0, 2.0].into_iter().collect(); + assert_eq!(e.median(), 1.0); + + // a single repeated value + let e: Empirical = [7.0; 5].into_iter().collect(); + assert_eq!(e.median(), 7.0); + + // one observation + let e: Empirical = [42.0].into_iter().collect(); + assert_eq!(e.median(), 42.0); + + // two observations straddle + let e: Empirical = [1.0, 4.0].into_iter().collect(); + assert_eq!(e.median(), 2.5); + + // negatives and duplicates together + let e: Empirical = [-5.0, -1.0, -1.0, 0.0, 3.0].into_iter().collect(); + assert_eq!(e.median(), -1.0); + } + + /// Agreement with the sorting reference across many shapes of data, + /// including heavy duplication, which is where the multiplicity handling in + /// the BTreeMap walk could go wrong. + #[test] + fn test_median_matches_sorting_reference() { + let cases: &[&[f64]] = &[ + &[1.0], + &[1.0, 2.0], + &[2.0, 1.0], + &[1.0, 2.0, 3.0], + &[1.0, 1.0, 2.0, 2.0], + &[1.0, 1.0, 1.0, 2.0, 2.0], + &[1.0, 2.0, 2.0, 2.0, 3.0], + &[5.0, 5.0, 5.0, 5.0], + &[1.0, 1.0, 1.0, 1.0, 9.0], + &[9.0, 1.0, 1.0, 1.0, 1.0], + &[-3.0, -2.0, -1.0, 0.0, 1.0, 2.0], + &[0.0, 0.0, 0.0, 1.0], + &[0.0, 1.0, 1.0, 1.0], + &[1e300, -1e300], + &[1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0], + ]; + for case in cases { + let e: Empirical = case.iter().copied().collect(); + let want = median_by_sorting(case); + let got = e.median(); + assert_eq!(got, want, "median of {case:?} was {got}, expected {want}"); + } + + // Longer runs built deterministically, with the modulus chosen to force + // many repeats. + for len in 1..40usize { + for modulus in [1u64, 2, 3, 7] { + let data: Vec = (0..len as u64) + .map(|i| ((i * 37 + 11) % modulus.max(1)) as f64) + .collect(); + let e: Empirical = data.iter().copied().collect(); + assert_eq!( + e.median(), + median_by_sorting(&data), + "len {len}, modulus {modulus}, data {data:?}" + ); + } + } + } + + /// The median must fall between the extremes, and coincide with them when + /// the data is constant. + #[test] + fn test_median_within_bounds() { + let e: Empirical = [1.0, 5.0, 2.0, 9.0, 3.0].into_iter().collect(); + assert!(e.median() >= e.min() && e.median() <= e.max()); + + let e: Empirical = [4.0; 3].into_iter().collect(); + assert_eq!(e.median(), e.min()); + assert_eq!(e.median(), e.max()); + } + + #[test] + #[should_panic(expected = "Cannot compute the median of zero samples")] + fn test_median_of_empty_panics() { + Empirical::new().unwrap().median(); + } + #[test] fn test_add_nan() { let mut empirical = Empirical::new().unwrap(); diff --git a/src/distribution/multivariate_normal.rs b/src/distribution/multivariate_normal.rs index e35a0fbd..312b8fb8 100644 --- a/src/distribution/multivariate_normal.rs +++ b/src/distribution/multivariate_normal.rs @@ -1,5 +1,5 @@ use crate::distribution::Continuous; -use crate::statistics::{Max, MeanN, Min, Mode, VarianceN}; +use crate::statistics::{Max, MeanN, Median, Min, Mode, VarianceN}; use core::f64; use core::f64::consts::{E, PI}; use nalgebra::{Cholesky, Const, DMatrix, DVector, Dim, DimMin, Dyn, OMatrix, OVector}; @@ -475,6 +475,36 @@ where } } +impl Median> for MultivariateNormal +where + D: Dim, + nalgebra::DefaultAllocator: + nalgebra::allocator::Allocator + nalgebra::allocator::Allocator, +{ + /// Returns the median of the multivariate normal distribution + /// + /// # Formula + /// + /// ```text + /// μ + /// ``` + /// + /// where `μ` is the mean + /// + /// # Remarks + /// + /// A multivariate distribution has several competing notions of median -- + /// the vector of marginal medians, the geometric (spatial) median, the + /// halfspace median -- which in general disagree. For a distribution + /// symmetric about a point they all coincide there, so for the multivariate + /// normal every one of them equals `μ`, as do the mean and the mode. That + /// agreement is what makes this unambiguous; it does not hold for + /// `Multinomial` or `Dirichlet`, which is why neither implements this trait. + fn median(&self) -> OVector { + self.mu.clone() + } +} + impl Mode> for MultivariateNormal where D: Dim, @@ -536,7 +566,7 @@ mod tests { use crate::{ distribution::{Continuous, MultivariateNormal}, - statistics::{Max, MeanN, Min, Mode, VarianceN}, + statistics::{Max, MeanN, Median, Min, Mode, VarianceN}, }; use super::MultivariateNormalError; @@ -672,6 +702,22 @@ mod tests { ); } + #[test] + fn test_median() { + let median = |x: MultivariateNormal<_>| x.median(); + test_case(vector![0., 0.], matrix![1., 0.; 0., 1.], vector![0., 0.], median); + test_case(vector![-3., 5.], matrix![2., 0.5; 0.5, 4.], vector![-3., 5.], median); + } + + /// Every notion of centre coincides for a symmetric distribution, which is + /// what makes `median` unambiguous here. + #[test] + fn test_median_mean_and_mode_coincide() { + let d = MultivariateNormal::new(vec![1.5, -2.0], vec![3.0, 0.7, 0.7, 2.0]).unwrap(); + assert_eq!(d.median(), d.mode()); + assert_eq!(d.median(), d.mean().unwrap()); + } + #[test] fn test_mode() { let mode = |x: MultivariateNormal<_>| x.mode(); diff --git a/src/distribution/multivariate_students_t.rs b/src/distribution/multivariate_students_t.rs index 2ad69880..03a1495d 100644 --- a/src/distribution/multivariate_students_t.rs +++ b/src/distribution/multivariate_students_t.rs @@ -1,6 +1,6 @@ use crate::distribution::Continuous; use crate::function::gamma; -use crate::statistics::{Max, MeanN, Min, Mode, VarianceN}; +use crate::statistics::{Max, MeanN, Median, Min, Mode, VarianceN}; use core::f64::consts::PI; use nalgebra::{Cholesky, Const, DMatrix, Dim, DimMin, Dyn, OMatrix, OVector}; @@ -305,6 +305,36 @@ where } } +impl Median> for MultivariateStudent +where + D: Dim, + nalgebra::DefaultAllocator: + nalgebra::allocator::Allocator + nalgebra::allocator::Allocator, +{ + /// Returns the median of the multivariate student's t-distribution + /// + /// # Formula + /// + /// ```text + /// μ + /// ``` + /// + /// where `μ` is the location + /// + /// # Remarks + /// + /// The distribution is symmetric about its location for every `ν`, so the + /// marginal, geometric and halfspace medians all coincide there; see + /// [`MultivariateNormal::median`](crate::distribution::MultivariateNormal::median). + /// + /// Note this is defined for `ν <= 1`, where + /// [`mean`](crate::statistics::MeanN::mean) is `None` because the first + /// moment does not converge. A median requires no moments to exist. + fn median(&self) -> OVector { + self.location.clone() + } +} + impl Mode> for MultivariateStudent where D: Dim, @@ -400,7 +430,7 @@ mod tests { use crate::{ distribution::{Continuous, MultivariateStudent, MultivariateNormal}, - statistics::{Max, MeanN, Min, Mode, VarianceN}, + statistics::{Max, MeanN, Median, Min, Mode, VarianceN}, }; use super::MultivariateStudentError; @@ -523,6 +553,26 @@ mod tests { test_case(vec![0., 0.], vec![1., 0., 0., 1.], 2., None, variance); } + #[test] + fn test_median() { + let median = |x: MultivariateStudent| x.median(); + test_case(vec![0., 0.], vec![1., 0., 0., 1.], 1., dvec![0., 0.], median); + test_case(vec![-1., 2.], vec![2., 0., 0., 3.], 5., dvec![-1., 2.], median); + } + + /// The median exists for every `freedom`, including the values at or below + /// one where the mean does not converge and `mean()` is `None`. + #[test] + fn test_median_defined_where_mean_is_not() { + for freedom in [0.5, 1.0] { + let d = MultivariateStudent::new(vec![3., -4.], vec![1., 0., 0., 1.], freedom).unwrap(); + assert!(d.mean().is_none(), "premise: mean undefined for v = {freedom}"); + assert_eq!(d.median(), dvec![3., -4.]); + // and it still agrees with the mode, by symmetry + assert_eq!(d.median(), d.mode()); + } + } + #[test] fn test_mode() { let mode = |x: MultivariateStudent| x.mode();