From 6500fd6556ee44bf6c2deee421d0cf26b4630655 Mon Sep 17 00:00:00 2001 From: Felix Agene Date: Sun, 26 Jul 2026 18:09:43 -0500 Subject: [PATCH 1/4] 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 810777fb..dbf8ed4a 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 alloc::{vec, vec::Vec}; -use nalgebra::{Dim, Dyn, OMatrix, OVector}; +use nalgebra::{Const, Dim, Dyn, OMatrix, OVector}; #[cfg(not(feature = "std"))] use num_traits::Float as _; @@ -217,6 +217,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, @@ -443,6 +482,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 5eaab233..0244389d 100644 --- a/src/distribution/multinomial.rs +++ b/src/distribution/multinomial.rs @@ -2,7 +2,7 @@ use crate::distribution::Discrete; use crate::function::factorial; use crate::statistics::*; use alloc::vec::Vec; -use nalgebra::{Dim, Dyn, OMatrix, OVector}; +use nalgebra::{Const, Dim, Dyn, OMatrix, OVector}; #[cfg(not(feature = "std"))] use num_traits::Float as _; @@ -267,6 +267,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, @@ -412,7 +462,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}; @@ -481,6 +531,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 91a016515aff076d7afead9bed79a448cd8a9423 Mon Sep 17 00:00:00 2001 From: Felix Agene Date: Sun, 26 Jul 2026 18:27:29 -0500 Subject: [PATCH 2/4] 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 38ae1517..c03e9828 100644 --- a/src/distribution/beta.rs +++ b/src/distribution/beta.rs @@ -320,6 +320,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 c4de36f0..851d29c3 100644 --- a/src/distribution/categorical.rs +++ b/src/distribution/categorical.rs @@ -312,6 +312,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` @@ -396,6 +421,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 208047db..a4dbf91a 100644 --- a/src/distribution/chi.rs +++ b/src/distribution/chi.rs @@ -300,6 +300,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 dbf8ed4a..cf0230b2 100644 --- a/src/distribution/dirichlet.rs +++ b/src/distribution/dirichlet.rs @@ -256,6 +256,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, @@ -323,13 +361,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 /// @@ -357,12 +402,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 /// @@ -386,25 +436,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) } } @@ -482,6 +530,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]); @@ -508,10 +626,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(); @@ -642,18 +762,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] @@ -663,20 +795,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 334e19b2..32231bde 100644 --- a/src/distribution/fisher_snedecor.rs +++ b/src/distribution/fisher_snedecor.rs @@ -210,6 +210,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 6df26e38..c3678c5b 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 ad45edb0..bd9672e3 100644 --- a/src/distribution/hypergeometric.rs +++ b/src/distribution/hypergeometric.rs @@ -368,6 +368,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 3630b072..0b3295f1 100644 --- a/src/distribution/inverse_gamma.rs +++ b/src/distribution/inverse_gamma.rs @@ -313,6 +313,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 fff79190..6e478f79 100644 --- a/src/distribution/negative_binomial.rs +++ b/src/distribution/negative_binomial.rs @@ -248,6 +248,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 5d54f41aba7140fb3b79b5fd63abf4e0cbcc807b Mon Sep 17 00:00:00 2001 From: Felix Agene Date: Sun, 26 Jul 2026 18:39:37 -0500 Subject: [PATCH 3/4] 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 6075a454..00b633bb 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 b0434eef..0e2f7d40 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 alloc::vec::Vec; use core::f64::consts::{E, PI}; use nalgebra::{Cholesky, Const, DMatrix, DVector, Dim, DimMin, Dyn, OMatrix, OVector}; @@ -477,6 +477,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, @@ -538,7 +568,7 @@ mod tests { use crate::{ distribution::{Continuous, MultivariateNormal}, - statistics::{Max, MeanN, Min, Mode, VarianceN}, + statistics::{Max, MeanN, Median, Min, Mode, VarianceN}, }; use super::MultivariateNormalError; @@ -674,6 +704,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 2423b8c2..647ab556 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 alloc::vec::Vec; use core::f64::consts::PI; use nalgebra::{Cholesky, Const, DMatrix, Dim, DimMin, Dyn, OMatrix, OVector}; @@ -308,6 +308,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, @@ -403,7 +433,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; @@ -526,6 +556,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(); From 3da838f6351002c04066a64a5b42b762a0943f7a Mon Sep 17 00:00:00 2001 From: "autofix-ci[bot]" <114827586+autofix-ci[bot]@users.noreply.github.com> Date: Wed, 26 Aug 2026 17:34:31 +0000 Subject: [PATCH 4/4] [autofix.ci] apply automated fixes --- src/distribution/empirical.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/distribution/empirical.rs b/src/distribution/empirical.rs index 00b633bb..6feb8a9d 100644 --- a/src/distribution/empirical.rs +++ b/src/distribution/empirical.rs @@ -258,7 +258,7 @@ impl Median for Empirical { // 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 need_two = n.is_multiple_of(2); let mut seen = 0; let mut lower = None;