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 5f1697b4..49ec338d 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,83 @@ 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 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, @@ -282,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 /// @@ -316,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 /// @@ -345,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) } } @@ -441,6 +528,115 @@ 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]); + 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. 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(); + 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(); @@ -564,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] @@ -585,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/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/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/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(); 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(); 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()); +}