From c8ed456c8697897e7bc44d6f4fc4c92d5b09a01e Mon Sep 17 00:00:00 2001
From: Orion Yeung <11580988+orionyeung001@users.noreply.github.com>
Date: Thu, 24 Oct 2024 12:29:00 -0500
Subject: [PATCH] fix: variance implementation for slices

bug introduced in a3f78fe554f18715ccc4ea906ff537a32dff1ad0
---
 src/statistics/slice_statistics.rs | 23 +++++++++++++++++------
 1 file changed, 17 insertions(+), 6 deletions(-)

diff --git a/src/statistics/slice_statistics.rs b/src/statistics/slice_statistics.rs
index b90b3cdf..74919f90 100644
--- a/src/statistics/slice_statistics.rs
+++ b/src/statistics/slice_statistics.rs
@@ -381,18 +381,29 @@ impl<D: AsMut<[f64]> + AsRef<[f64]>> StandardizedMoment<f64> for Data<D> {
         // TODO: rewrite iterator statistics and rely on them instead
         let mut iter = self.0.as_ref().iter().enumerate();
 
-        let mut sum = *iter.next()?.1;
-        let mut variance = 0.0;
+        // init single sample (mean=sample and M_2=0)
+        let mut n: f64 = 1.;
+        let mut mu = *iter.next()?.1;
+        let mut m2 = 0.0;
 
         for (i, &x) in iter {
             if x.is_nan() {
+                // fail early
                 return Some(f64::NAN);
             }
-            let diff = (i + 1) as f64 * x - sum;
-            sum += x;
-            variance += diff * diff / (i * (i + 1)) as f64
+            n = (i + 1) as f64;
+            let delta = x - mu;
+            mu += delta / n;
+            m2 += delta * (x - mu);
+        }
+
+        if n == 1. {
+            // no sample variance for a single sample
+            None
+        } else {
+            // normalize by bessel
+            Some(m2 / (n - 1.))
         }
-        Some(variance / (self.0.as_ref().len() - 1) as f64)
     }
 
     fn skewness(&self) -> Self::Skew {