Merge pull request #486 from LukeMathWalker/std_axis

Author: jturner314

src/numeric/impl_numeric.rs

@@ -181,6 +181,53 @@ impl<A, S, D> ArrayBase<S, D>
         }
     }
 
+    /// Return standard deviation along `axis`.
+    ///
+    /// The standard deviation is computed from the variance using
+    /// the [Welford one-pass algorithm](https://www.jstor.org/stable/1266577).
+    ///
+    /// The parameter `ddof` specifies the "delta degrees of freedom". For
+    /// example, to calculate the population standard deviation, use `ddof = 0`,
+    /// or to calculate the sample standard deviation, use `ddof = 1`.
+    ///
+    /// The standard deviation is defined as:
+    ///
+    /// ```text
+    ///                    1       n
+    /// stddev = sqrt ( ――――――――   ∑ (xᵢ - x̅)² )
+    ///                 n - ddof  i=1
+    /// ```
+    ///
+    /// where
+    ///
+    /// ```text
+    ///     1   n
+    /// x̅ = ―   ∑ xᵢ
+    ///     n  i=1
+    /// ```
+    ///
+    /// **Panics** if `ddof` is greater than or equal to the length of the
+    /// axis, if `axis` is out of bounds, or if the length of the axis is zero.
+    ///
+    /// # Example
+    ///
+    /// ```
+    /// use ndarray::{aview1, arr2, Axis};
+    ///
+    /// let a = arr2(&[[1., 2.],
+    ///                [3., 4.],
+    ///                [5., 6.]]);
+    /// let stddev = a.std_axis(Axis(0), 1.);
+    /// assert_eq!(stddev, aview1(&[2., 2.]));
+    /// ```
+    pub fn std_axis(&self, axis: Axis, ddof: A) -> Array<A, D::Smaller>
+    where
+        A: Float,
+        D: RemoveAxis,
+    {
+        self.var_axis(axis, ddof).mapv_into(|x| x.sqrt())
+    }
+
     /// Return `true` if the arrays' elementwise differences are all within
     /// the given absolute tolerance, `false` otherwise.
     ///

tests/array.rs

@@ -745,6 +745,53 @@ fn var_axis() {
     assert!(d.var_axis(Axis(0), 0.).all_close(&aview0(&1.8875), 1e-12));
 }
 
+#[test]
+fn std_axis() {
+    let a = array![
+       [
+            [ 0.22935481,  0.08030619,  0.60827517,  0.73684379],
+            [ 0.90339851,  0.82859436,  0.64020362,  0.2774583 ],
+            [ 0.44485313,  0.63316367,  0.11005111,  0.08656246]
+       ],
+       [
+           [ 0.28924665,  0.44082454,  0.59837736,  0.41014531],
+            [ 0.08382316,  0.43259439,  0.1428889 ,  0.44830176],
+            [ 0.51529756,  0.70111616,  0.20799415,  0.91851457]
+       ],
+    ];
+    assert!(a.std_axis(Axis(0), 1.5).all_close(
+        &aview2(&[
+            [ 0.05989184,  0.36051836,  0.00989781,  0.32669847],
+            [ 0.81957535,  0.39599997,  0.49731472,  0.17084346],
+            [ 0.07044443,  0.06795249,  0.09794304,  0.83195211],
+        ]),
+        1e-4,
+    ));
+    assert!(a.std_axis(Axis(1), 1.7).all_close(
+        &aview2(&[
+            [ 0.42698655,  0.48139215,  0.36874991,  0.41458724],
+            [ 0.26769097,  0.18941435,  0.30555015,  0.35118674],
+        ]),
+        1e-8,
+    ));
+    assert!(a.std_axis(Axis(2), 2.3).all_close(
+        &aview2(&[
+            [ 0.41117907,  0.37130425,  0.35332388],
+            [ 0.16905862,  0.25304841,  0.39978276],
+        ]),
+        1e-8,
+    ));
+
+    let b = array![[100000., 1., 0.01]];
+    assert!(b.std_axis(Axis(0), 0.).all_close(&aview1(&[0., 0., 0.]), 1e-12));
+    assert!(
+        b.std_axis(Axis(1), 0.).all_close(&aview1(&[47140.214021552769]), 1e-6),
+    );
+
+    let c = array![[], []];
+    assert_eq!(c.std_axis(Axis(0), 0.), aview1(&[]));
+}
+
 #[test]
 #[should_panic]
 fn var_axis_bad_dof() {
@@ -759,6 +806,20 @@ fn var_axis_empty_axis() {
     a.var_axis(Axis(1), 0.);
 }
 
+#[test]
+#[should_panic]
+fn std_axis_bad_dof() {
+    let a = array![1., 2., 3.];
+    a.std_axis(Axis(0), 4.);
+}
+
+#[test]
+#[should_panic]
+fn std_axis_empty_axis() {
+    let a = array![[], []];
+    a.std_axis(Axis(1), 0.);
+}
+
 #[test]
 fn iter_size_hint()
 {