From 46e4beb83f345500c696040f427715164aa02597 Mon Sep 17 00:00:00 2001
From: oberien <jaro.fietz@gmx.de>
Date: Wed, 1 May 2019 12:38:51 +0200
Subject: [PATCH] Add standard deviation

---
 src/summary_statistics/means.rs | 19 +++++++++++++++++++
 src/summary_statistics/mod.rs   | 18 ++++++++++++++++++
 2 files changed, 37 insertions(+)

diff --git a/src/summary_statistics/means.rs b/src/summary_statistics/means.rs
index 89d4df9d..c42a36bf 100644
--- a/src/summary_statistics/means.rs
+++ b/src/summary_statistics/means.rs
@@ -54,6 +54,13 @@ where
         Ok(central_moments[3] / central_moments[2].sqrt().powi(3))
     }
 
+    fn standard_deviation(&self) -> Result<A, EmptyInput>
+    where
+        A: Float + FromPrimitive,
+    {
+        Ok(self.central_moment(2)?.sqrt())
+    }
+
     fn central_moment(&self, order: u16) -> Result<A, EmptyInput>
     where
         A: Float + FromPrimitive,
@@ -277,6 +284,18 @@ mod tests {
         assert_eq!(a.central_moment(1).unwrap(), 0.);
     }
 
+    #[test]
+    fn test_standard_deviation() {
+        let a: Array1<f64> = array![4.0, 9.0, 11.0, 12.0, 17.0, 5.0, 8.0, 12.0, 14.0];
+        // Computed using numpy.std
+        let expected_std = 3.9377878103709665;
+        assert_abs_diff_eq!(
+            a.standard_deviation().unwrap(),
+            expected_std,
+            epsilon = 1e-12
+        );
+    }
+
     #[test]
     fn test_central_moments() {
         let a: Array1<f64> = array![
diff --git a/src/summary_statistics/mod.rs b/src/summary_statistics/mod.rs
index 8351ff82..5a9d9426 100644
--- a/src/summary_statistics/mod.rs
+++ b/src/summary_statistics/mod.rs
@@ -101,6 +101,24 @@ where
     where
         A: Float + FromPrimitive;
 
+    /// Returns the [standard deviation] σ of all elements in the array:
+    /// ```text
+    ///       _____________
+    ///      / 1  n
+    /// σ = /  ―  ∑ (xᵢ-x̅)²
+    ///    √   n i=1
+    /// ```
+    ///
+    /// If the array is empty, `Err(EmptyInput)` is returned.
+    ///
+    /// **Panics** if `A::from_usize()` fails to convert the number of elements
+    /// in the array or if `order` overflows `i32`.
+    ///
+    /// [standard deviation]: https://en.wikipedia.org/wiki/Standard_deviation
+    fn standard_deviation(&self) -> Result<A, EmptyInput>
+    where
+        A: Float + FromPrimitive;
+
     /// Returns the *p*-th [central moment] of all elements in the array, μₚ:
     ///
     /// ```text