feat: Adding Mode and Moments to summary statistics

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.);

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.);

src/summary_statistics/means.rs

@@ -104,6 +104,89 @@ where
             .ok_or(EmptyInput)
     }
 
+    fn mode(&self) -> Result<A, EmptyInput>
+    where
+        A: Clone + PartialEq,
+    {
+        self.modes()?.into_iter().next().ok_or(EmptyInput)
+    }
+
+    fn modes(&self) -> Result<Vec<A>, EmptyInput>
+    where
+        A: Clone + PartialEq,
+    {
+        if self.is_empty() {
+            return Err(EmptyInput);
+        }
+        Ok(modes(self.iter().cloned()))
+    }
+
+    fn mode_axis(&self, axis: Axis) -> Result<Array<A, D::Smaller>, 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<A, EmptyInput>
+    where
+        A: Float + FromPrimitive,
+    {
+        self.raw_moments(order)
+            .map(|moments| moments[usize::from(order)])
+    }
+
+    fn raw_moments(&self, order: u16) -> Result<Vec<A>, EmptyInput>
+    where
+        A: Float + FromPrimitive,
+    {
+        if self.is_empty() {
+            return Err(EmptyInput);
+        }
+        Ok(moments(self.view(), order))
+    }
+
+    fn standardized_moment(&self, order: u16) -> Result<A, EmptyInput>
+    where
+        A: Float + FromPrimitive,
+    {
+        self.standardized_moments(order)
+            .map(|moments| moments[usize::from(order)])
+    }
+
+    fn standardized_moments(&self, order: u16) -> Result<Vec<A>, 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<A, MultiInputError>
     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<A, I>(values: I) -> Vec<A>
+where
+    A: Clone + PartialEq,
+    I: IntoIterator<Item = A>,
+{
+    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:
 ///

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<A, EmptyInput>
+    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<Vec<A>, 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<Array<A, D::Smaller>, 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<A, EmptyInput>
+    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<Vec<A>, 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<A, EmptyInput>
+    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<Vec<A>, 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.

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<f64> = 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);
+}