diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 6fc4f6cf4..a3858f829 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -13,6 +13,8 @@ use crate::imp_prelude::*;
 use crate::itertools::enumerate;
 use crate::numeric_util;
 
+use crate::{FoldWhile, Zip};
+
 /// # Numerical Methods for Arrays
 impl<A, S, D> ArrayBase<S, D>
 where
@@ -111,6 +113,114 @@ where
         sum
     }
 
+    /// Return variance of elements in the array.
+    ///
+    /// The variance is computed 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 variance, use `ddof = 0`, or to
+    /// calculate the sample variance, use `ddof = 1`.
+    ///
+    /// The variance is defined as:
+    ///
+    /// ```text
+    ///               1       n
+    /// variance = ――――――――   ∑ (xᵢ - x̅)²
+    ///            n - ddof  i=1
+    /// ```
+    ///
+    /// where
+    ///
+    /// ```text
+    ///     1   n
+    /// x̅ = ―   ∑ xᵢ
+    ///     n  i=1
+    /// ```
+    ///
+    /// and `n` is the length of the array.
+    ///
+    /// **Panics** if `ddof` is less than zero or greater than `n`
+    ///
+    /// # Example
+    ///
+    /// ```
+    /// use ndarray::array;
+    /// use approx::assert_abs_diff_eq;
+    ///
+    /// let a = array![1., -4.32, 1.14, 0.32];
+    /// let var = a.var(1.);
+    /// assert_abs_diff_eq!(var, 6.7331, epsilon = 1e-4);
+    /// ```
+    pub fn var(&self, ddof: A) -> A
+    where
+        A: Float + FromPrimitive,
+    {
+        let zero = A::from_usize(0).expect("Converting 0 to `A` must not fail.");
+        let n = A::from_usize(self.len()).expect("Converting length to `A` must not fail.");
+        assert!(
+            !(ddof < zero || ddof > n),
+            "`ddof` must not be less than zero or greater than the length of \
+             the axis",
+        );
+        let dof = n - ddof;
+        let mut mean = A::zero();
+        let mut sum_sq = A::zero();
+        for (i, &x) in self.into_iter().enumerate() {
+            let count = A::from_usize(i + 1).expect("Converting index to `A` must not fail.");
+            let delta = x - mean;
+            mean = mean + delta / count;
+            sum_sq = (x - mean).mul_add(delta, sum_sq);
+        }
+        sum_sq / dof
+    }
+
+    /// Return standard deviation of elements in the array.
+    ///
+    /// 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
+    /// ```
+    ///
+    /// and `n` is the length of the array.
+    ///
+    /// **Panics** if `ddof` is less than zero or greater than `n`
+    ///
+    /// # Example
+    ///
+    /// ```
+    /// use ndarray::array;
+    /// use approx::assert_abs_diff_eq;
+    ///
+    /// let a = array![1., -4.32, 1.14, 0.32];
+    /// let stddev = a.std(1.);
+    /// assert_abs_diff_eq!(stddev, 2.59483, epsilon = 1e-4);
+    /// ```
+    pub fn std(&self, ddof: A) -> A
+    where
+        A: Float + FromPrimitive,
+    {
+        self.var(ddof).sqrt()
+    }
+
     /// Return sum along `axis`.
     ///
     /// ```
diff --git a/src/private.rs b/src/private.rs
index ea13164e4..552b31497 100644
--- a/src/private.rs
+++ b/src/private.rs
@@ -21,5 +21,5 @@ macro_rules! private_impl {
         fn __private__(&self) -> crate::private::PrivateMarker {
             crate::private::PrivateMarker
         }
-    }
+    };
 }
diff --git a/tests/numeric.rs b/tests/numeric.rs
index 7c6f1441e..13f8433c5 100644
--- a/tests/numeric.rs
+++ b/tests/numeric.rs
@@ -64,6 +64,72 @@ fn sum_mean_empty() {
     assert_eq!(a, None);
 }
 
+#[test]
+fn var() {
+    let a = array![1., -4.32, 1.14, 0.32];
+    assert_abs_diff_eq!(a.var(0.), 5.049875, epsilon = 1e-8);
+}
+
+#[test]
+#[should_panic]
+fn var_negative_ddof() {
+    let a = array![1., 2., 3.];
+    a.var(-1.);
+}
+
+#[test]
+#[should_panic]
+fn var_too_large_ddof() {
+    let a = array![1., 2., 3.];
+    a.var(4.);
+}
+
+#[test]
+fn var_nan_ddof() {
+    let a = Array2::<f64>::zeros((2, 3));
+    let v = a.var(::std::f64::NAN);
+    assert!(v.is_nan());
+}
+
+#[test]
+fn var_empty_arr() {
+    let a: Array1<f64> = array![];
+    assert!(a.var(0.0).is_nan());
+}
+
+#[test]
+fn std() {
+    let a = array![1., -4.32, 1.14, 0.32];
+    assert_abs_diff_eq!(a.std(0.), 2.24719, epsilon = 1e-5);
+}
+
+#[test]
+#[should_panic]
+fn std_negative_ddof() {
+    let a = array![1., 2., 3.];
+    a.std(-1.);
+}
+
+#[test]
+#[should_panic]
+fn std_too_large_ddof() {
+    let a = array![1., 2., 3.];
+    a.std(4.);
+}
+
+#[test]
+fn std_nan_ddof() {
+    let a = Array2::<f64>::zeros((2, 3));
+    let v = a.std(::std::f64::NAN);
+    assert!(v.is_nan());
+}
+
+#[test]
+fn std_empty_arr() {
+    let a: Array1<f64> = array![];
+    assert!(a.std(0.0).is_nan());
+}
+
 #[test]
 #[cfg(feature = "approx")]
 fn var_axis() {