diff --git a/README.md b/README.md index dc53029..0128785 100644 --- a/README.md +++ b/README.md @@ -9,7 +9,7 @@ This crate provides statistical methods for [`ndarray`]'s `ArrayRef` type. Currently available routines include: - order statistics (minimum, maximum, median, quantiles, etc.); -- summary statistics (mean, skewness, kurtosis, central moments, etc.) +- summary statistics (mean, mode, raw/central/standardized moments, skewness, kurtosis, etc.) - partitioning; - correlation analysis (covariance, pearson correlation); - measures from information theory (entropy, KL divergence, etc.); diff --git a/src/lib.rs b/src/lib.rs index ef4c5c4..f824e6d 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -3,7 +3,7 @@ //! //! Currently available routines include: //! - [order statistics] (minimum, maximum, median, quantiles, etc.); -//! - [summary statistics] (mean, skewness, kurtosis, central moments, etc.) +//! - [summary statistics] (mean, mode, raw/central/standardized moments, skewness, kurtosis, etc.) //! - [partitioning]; //! - [correlation analysis] (covariance, pearson correlation); //! - [measures from information theory] (entropy, KL divergence, etc.); diff --git a/src/summary_statistics/means.rs b/src/summary_statistics/means.rs index 92cd101..a2e8ab0 100644 --- a/src/summary_statistics/means.rs +++ b/src/summary_statistics/means.rs @@ -104,6 +104,89 @@ where .ok_or(EmptyInput) } + fn mode(&self) -> Result + where + A: Clone + PartialEq, + { + self.modes()?.into_iter().next().ok_or(EmptyInput) + } + + fn modes(&self) -> Result, EmptyInput> + where + A: Clone + PartialEq, + { + if self.is_empty() { + return Err(EmptyInput); + } + Ok(modes(self.iter().cloned())) + } + + fn mode_axis(&self, axis: Axis) -> Result, EmptyInput> + where + A: Clone + PartialEq, + D: RemoveAxis, + { + if self.is_empty() { + return Err(EmptyInput); + } + + Ok(self.map_axis(axis, |lane| { + modes(lane.iter().cloned()) + .into_iter() + .next() + .expect("non-empty lanes must have a mode") + })) + } + + fn raw_moment(&self, order: u16) -> Result + where + A: Float + FromPrimitive, + { + self.raw_moments(order) + .map(|moments| moments[usize::from(order)]) + } + + fn raw_moments(&self, order: u16) -> Result, EmptyInput> + where + A: Float + FromPrimitive, + { + if self.is_empty() { + return Err(EmptyInput); + } + Ok(moments(self.view(), order)) + } + + fn standardized_moment(&self, order: u16) -> Result + where + A: Float + FromPrimitive, + { + self.standardized_moments(order) + .map(|moments| moments[usize::from(order)]) + } + + fn standardized_moments(&self, order: u16) -> Result, EmptyInput> + where + A: Float + FromPrimitive, + { + let central_moments = self.central_moments(order)?; + if order < 2 { + return Ok(central_moments); + } + + let standard_deviation = central_moments[2].sqrt(); + Ok(central_moments + .into_iter() + .enumerate() + .map(|(order, moment)| { + if order < 2 { + moment + } else { + moment / standard_deviation.powi(order as i32) + } + }) + .collect()) + } + fn weighted_var(&self, weights: &Self, ddof: A) -> Result where A: AddAssign + Float + FromPrimitive, @@ -266,6 +349,35 @@ where Ok(s / (weight_sum - ddof)) } +/// Returns all values with the greatest frequency, preserving first +/// occurrence order. This deliberately uses `PartialEq` instead of `Hash` so +/// that ordinary floating-point arrays can use the mode API. +fn modes(values: I) -> Vec +where + A: Clone + PartialEq, + I: IntoIterator, +{ + let mut counts: Vec<(A, usize)> = Vec::new(); + for value in values { + if let Some(index) = counts.iter().position(|(candidate, _)| candidate == &value) { + counts[index].1 += 1; + } else { + counts.push((value, 1)); + } + } + + let greatest_frequency = counts + .iter() + .map(|(_, frequency)| *frequency) + .max() + .unwrap_or(0); + counts + .into_iter() + .filter(|(_, frequency)| *frequency == greatest_frequency) + .map(|(value, _)| value) + .collect() +} + /// Returns a vector containing all moments of the array elements up to /// *order*, where the *p*-th moment is defined as: /// diff --git a/src/summary_statistics/mod.rs b/src/summary_statistics/mod.rs index e31dde1..239440e 100644 --- a/src/summary_statistics/mod.rs +++ b/src/summary_statistics/mod.rs @@ -155,6 +155,90 @@ where where A: Float + FromPrimitive; + /// Returns the mode of all elements in the array. + /// + /// If several values have the same greatest frequency, the value whose + /// first occurrence is earliest in the array is returned. Values are + /// compared with `PartialEq`, so this method also supports floating-point + /// arrays without requiring `Eq` or `Hash`. In particular, NaN values do + /// not compare equal to one another. + /// + /// If the array is empty, `Err(EmptyInput)` is returned. + fn mode(&self) -> Result + where + A: Clone + PartialEq; + + /// Returns all modes of the array, ordered by their first occurrence. + /// + /// If several values have the same greatest frequency, all of them are + /// returned. Values are compared with `PartialEq`; see [`mode`] for the + /// implications for floating-point NaN values. + /// + /// If the array is empty, `Err(EmptyInput)` is returned. + /// + /// [`mode`]: #tymethod.mode + fn modes(&self) -> Result, EmptyInput> + where + A: Clone + PartialEq; + + /// Returns the mode along `axis` for every one-dimensional lane. + /// + /// Ties are resolved by first occurrence within each lane. If the array + /// is empty, `Err(EmptyInput)` is returned. As with other axis methods, + /// this method panics if `axis` is out of bounds. + fn mode_axis(&self, axis: Axis) -> Result, EmptyInput> + where + A: Clone + PartialEq, + D: RemoveAxis; + + /// Returns the *p*-th raw moment of all elements in the array: + /// + /// ```text + /// 1 n + /// mₚ = ─ ∑ xᵢᵖ + /// n i=1 + /// ``` + /// + /// The zeroth raw moment is one and the first raw moment is the arithmetic + /// mean. If the array is empty, `Err(EmptyInput)` is returned. + fn raw_moment(&self, order: u16) -> Result + where + A: Float + FromPrimitive; + + /// Returns all raw moments from order zero through `order`. + /// + /// The returned vector is indexed by moment order. If the array is empty, + /// `Err(EmptyInput)` is returned. + fn raw_moments(&self, order: u16) -> Result, EmptyInput> + where + A: Float + FromPrimitive; + + /// Returns the *p*-th standardized moment: + /// + /// ```text + /// μₚ + /// ηₚ = ──── + /// σᵖ + /// ``` + /// + /// The zeroth and first standardized moments are one and zero, + /// respectively. The third and fourth standardized moments are skewness + /// and Pearson's kurtosis. For a zero-variance input, moments of order two + /// or greater follow floating-point division semantics and are generally + /// NaN. + /// + /// If the array is empty, `Err(EmptyInput)` is returned. + fn standardized_moment(&self, order: u16) -> Result + where + A: Float + FromPrimitive; + + /// Returns all standardized moments from order zero through `order`. + /// + /// If the array is empty, `Err(EmptyInput)` is returned. + fn standardized_moments(&self, order: u16) -> Result, EmptyInput> + where + A: Float + FromPrimitive; + /// Return weighted variance of all elements in the array. /// /// The weighted variance is computed using the [`West, D. H. D.`] incremental algorithm. diff --git a/tests/summary_statistics.rs b/tests/summary_statistics.rs index 1062a09..af6256a 100644 --- a/tests/summary_statistics.rs +++ b/tests/summary_statistics.rs @@ -412,3 +412,61 @@ fn test_kurtosis_and_skewness() { assert_abs_diff_eq!(kurtosis, expected_kurtosis, epsilon = 1e-12); assert_abs_diff_eq!(skewness, expected_skewness, epsilon = 1e-8); } + +#[test] +fn test_mode_and_modes() { + let a = array![3.0, 1.0, 3.0, 2.0, 2.0, 4.0]; + + // 3 and 2 are tied; mode() resolves ties by first occurrence. + assert_eq!(a.mode().unwrap(), 3.0); + assert_eq!(a.modes().unwrap(), vec![3.0, 2.0]); +} + +#[test] +fn test_mode_axis() { + let a = array![[1, 2, 2, 3], [4, 5, 4, 5]]; + + assert_eq!(a.mode_axis(Axis(0)).unwrap(), array![1, 2, 2, 3]); + assert_eq!(a.mode_axis(Axis(1)).unwrap(), array![2, 4]); +} + +#[test] +fn test_mode_and_moments_with_empty_array() { + let a: Array1 = array![]; + + assert_eq!(a.mode(), Err(EmptyInput)); + assert_eq!(a.modes(), Err(EmptyInput)); + assert_eq!(a.mode_axis(Axis(0)), Err(EmptyInput)); + assert_eq!(a.raw_moment(2), Err(EmptyInput)); + assert_eq!(a.raw_moments(2), Err(EmptyInput)); + assert_eq!(a.standardized_moment(2), Err(EmptyInput)); + assert_eq!(a.standardized_moments(2), Err(EmptyInput)); +} + +#[test] +fn test_raw_and_standardized_moments() { + let a = array![1.0, 2.0, 2.0, 3.0]; + + assert_abs_diff_eq!(a.raw_moment(0).unwrap(), 1.0, epsilon = 1e-12); + assert_abs_diff_eq!(a.raw_moment(1).unwrap(), 2.0, epsilon = 1e-12); + assert_abs_diff_eq!(a.raw_moment(2).unwrap(), 4.5, epsilon = 1e-12); + assert_abs_diff_eq!(a.raw_moment(3).unwrap(), 11.0, epsilon = 1e-12); + for (actual, expected) in a + .raw_moments(4) + .unwrap() + .into_iter() + .zip([1.0, 2.0, 4.5, 11.0, 28.5]) + { + assert_abs_diff_eq!(actual, expected, epsilon = 1e-12); + } + + for (actual, expected) in a + .standardized_moments(4) + .unwrap() + .into_iter() + .zip([1.0, 0.0, 1.0, 0.0, 2.0]) + { + assert_abs_diff_eq!(actual, expected, epsilon = 1e-12); + } + assert_abs_diff_eq!(a.standardized_moment(3).unwrap(), 0.0, epsilon = 1e-12); +}