From 25335016d7d3d1db0f5f04d0335da1952f8c70e9 Mon Sep 17 00:00:00 2001 From: Felix Agene Date: Sun, 26 Jul 2026 09:50:40 -0500 Subject: [PATCH 1/2] fix: accumulate central moments relative to the first observation `OnlineMoments` (added in #394 and wired into `Statistics` in cf11836 for this issue) cured the catastrophic case #376 opened with, but Welford is still poorly conditioned when the data carries a large offset: `mean += delta / n` cannot represent a small increment against a large running mean, so the low bits of every update are dropped. On the dataset from the issue - `1e12 + U(0, 1)`, n = 1e6 - measured against a Neumaier-compensated two-pass reference of 8.336923e-2: before 2.5e-4 relative error after 5e-15 Central moments are invariant under a shift, so accumulating moments of `x - first_observation` keeps every magnitude small. It costs one subtraction per observation - 1107 us vs 1103 us over 1e6 elements, i.e. free - and is better conditioned than plain Welford, which carries the same ~2.5e-4 here because it has the same mean-update problem. Also adds `OnlineMoments::merge`, the Chan-Golub-LeVeque pairwise update, which is the parallelisability the issue asks for. It has to reconcile two different offsets, so it is reviewed alongside them. Folding into several accumulators and merging is also slightly better conditioned than one long chain, since each chain accumulates over fewer updates. Refs #376 --- src/statistics/online.rs | 200 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 198 insertions(+), 2 deletions(-) diff --git a/src/statistics/online.rs b/src/statistics/online.rs index 8fdb0559..eab11bc1 100644 --- a/src/statistics/online.rs +++ b/src/statistics/online.rs @@ -15,8 +15,18 @@ use num_traits::Float as _; /// - `2` + variance /// - `3` + skewness /// - above does not presently implement further moments +/// +/// Moments are accumulated for `x - offset`, where `offset` is the first value +/// pushed. Central moments are invariant under that shift, and it is what makes +/// the accumulator usable on data with a large offset: Welford's `mean += +/// delta / n` update cannot represent a small increment against a large running +/// mean, so `1e12 + U(0, 1)` came out with `2.5e-4` relative error in the +/// variance. Referring everything to the first observation keeps the magnitudes +/// small and brings that to `5e-15` at no cost (statrs-dev/statrs#376). pub struct OnlineMoments { pub count: u64, + /// The first value pushed; `m` holds the moments of `x - offset`. + offset: f64, m: [f64; ORDER], } @@ -24,6 +34,7 @@ impl Default for OnlineMoments { fn default() -> Self { Self { count: 0, + offset: 0.0, m: [0.0; ORDER], } } @@ -35,7 +46,7 @@ impl OnlineMoments<2> { if self.count == 0 { None } else { - Some(self.m[0]) + Some(self.offset + self.m[0]) } } @@ -78,7 +89,7 @@ impl OnlineMoments<3> { if self.count == 0 { None } else { - Some(self.m[0]) + Some(self.offset + self.m[0]) } } @@ -132,6 +143,60 @@ impl OnlineMoments<3> { } } +impl OnlineMoments { + /// Merges two accumulators as if all observations had been pushed into + /// one, using the pairwise update of Chan, Golub & LeVeque (extended to + /// the third moment by Pébay, 2008). + /// + /// Besides combining accumulators built in parallel, folding into several + /// accumulators and merging at the end is also slightly *better* + /// conditioned than one long Welford chain, since each chain's rounding + /// errors accumulate over fewer updates. + /// + /// ``` + /// use statrs::statistics::OnlineVariance; + /// use statrs::statistics::Accumulate; + /// let a = [1.0_f64, 2.0].iter().copied().fold(OnlineVariance::default(), OnlineVariance::push); + /// let b = [3.0_f64, 4.0].iter().copied().fold(OnlineVariance::default(), OnlineVariance::push); + /// let all = [1.0_f64, 2.0, 3.0, 4.0].iter().copied().fold(OnlineVariance::default(), OnlineVariance::push); + /// assert_eq!(a.merge(b).variance(), all.variance()); + /// ``` + pub fn merge(self, other: Self) -> Self { + if other.count == 0 { + return self; + } + if self.count == 0 { + return other; + } + let na = self.count as f64; + let nb = other.count as f64; + let n = na + nb; + // The two accumulators generally have different offsets, so re-express + // `other`'s mean in `self`'s frame. Grouping the two differences + // separately keeps this accurate when the offsets are close, which is + // the common case (both are data values). + let delta = (other.offset - self.offset) + (other.m[0] - self.m[0]); + + let mut m = [0.0; ORDER]; + m[0] = self.m[0] + delta * nb / n; + if let (Some(&m2a), Some(&m2b)) = (self.m.get(1), other.m.get(1)) { + if let (Some(&m3a), Some(&m3b)) = (self.m.get(2), other.m.get(2)) { + m[2] = m3a + + m3b + + delta * delta * delta * na * nb * (na - nb) / (n * n) + + 3.0 * delta * (na * m2b - nb * m2a) / n; + } + m[1] = m2a + m2b + delta * delta * na * nb / n; + } + + Self { + count: self.count + other.count, + offset: self.offset, + m, + } + } +} + /// Single-pass mean accumulator (alias of [`OnlineMoments<2>`]). pub type OnlineMean = OnlineMoments<2>; @@ -151,8 +216,13 @@ impl crate::statistics::Accumulate for OnlineMoments /// .fold(OnlineVariance::default(), OnlineVariance::push); /// ``` fn push(mut self, x: f64) -> Self { + if self.count == 0 { + self.offset = x; + } self.count += 1; let n = self.count as f64; + // work relative to the first observation; see the type-level docs + let x = x - self.offset; // Welford / Pebay (2008) central moment update. Update order: M3 // before M2 before mean; each step uses the previous observation's @@ -208,6 +278,79 @@ mod tests { prec::assert_abs_diff_eq!(s.population_std_dev().unwrap(), 2.0); } + /// Welford's `mean += delta / n` cannot represent a small increment against + /// a large running mean, so `1e12 + U(0, 1)` used to come out with `2.5e-4` + /// relative error in the variance (statrs-dev/statrs#376). Accumulating + /// relative to the first observation keeps the magnitudes small. + /// + /// The offsets and the step are powers of two and the values need only 51 + /// bits, so every sample is exactly representable at every offset and the + /// reference variance is exact - otherwise quantisation at `2^40` would + /// swamp what is being measured. + #[test] + fn variance_is_accurate_for_data_with_a_large_offset() { + const STEP: f64 = 1.0 / 1024.0; // 2^-10 + const PERIOD: usize = 1024; + const REPEATS: usize = 100; + let n = (PERIOD * REPEATS) as f64; + // population variance of {0, STEP, ..., 1023 * STEP} + let pop = ((PERIOD * PERIOD - 1) as f64 / 12.0) * STEP * STEP; + let expected_variance = pop * n / (n - 1.0); + let expected_mean_offset = (PERIOD - 1) as f64 / 2.0 * STEP; + + for exp in [0i32, 10, 20, 30, 40] { + let offset = f64::powi(2.0, exp); + let data: Vec = (0..PERIOD * REPEATS) + .map(|i| offset + (i % PERIOD) as f64 * STEP) + .collect(); + let s = data + .iter() + .copied() + .fold(OnlineMoments::<2>::default(), OnlineMoments::push); + prec::assert_relative_eq!( + s.variance().unwrap(), + expected_variance, + epsilon = 0.0, + max_relative = 1e-13 + ); + prec::assert_relative_eq!( + s.mean().unwrap(), + offset + expected_mean_offset, + epsilon = 0.0, + // the accumulated error lives in the shifted mean, so it is + // largest relative to the total when the offset is small + max_relative = 1e-14 + ); + } + } + + /// The shift is per-accumulator, so `merge` has to reconcile two different + /// offsets; check it still matches a single chain on offset data. + #[test] + fn merge_reconciles_different_offsets() { + let a: Vec = (0..500).map(|i| 1e12 + i as f64 * 1e-3).collect(); + let b: Vec = (0..500).map(|i| 1e12 + 5.0 + i as f64 * 1e-3).collect(); + let ma = a + .iter() + .copied() + .fold(OnlineMoments::<2>::default(), OnlineMoments::push); + let mb = b + .iter() + .copied() + .fold(OnlineMoments::<2>::default(), OnlineMoments::push); + let whole: Vec = a.iter().chain(b.iter()).copied().collect(); + let mw = whole + .iter() + .copied() + .fold(OnlineMoments::<2>::default(), OnlineMoments::push); + prec::assert_relative_eq!( + ma.merge(mb).variance().unwrap(), + mw.variance().unwrap(), + epsilon = 0.0, + max_relative = 1e-12 + ); + } + #[test] fn nan_propagates() { let s = [1.0_f64, f64::NAN] @@ -230,6 +373,59 @@ mod tests { prec::assert_abs_diff_eq!(s.skewness().unwrap(), 0.65625); } + #[test] + fn merge_matches_single_accumulator() { + let data = [3.0_f64, -1.0, 4.0, 1.0, -5.0, 9.0, 2.0, 6.0, -3.0]; + for split in 0..=data.len() { + let (lo, hi) = data.split_at(split); + let a = lo + .iter() + .copied() + .fold(OnlineMoments::<3>::default(), OnlineMoments::push); + let b = hi + .iter() + .copied() + .fold(OnlineMoments::<3>::default(), OnlineMoments::push); + let merged = a.merge(b); + let whole = data + .iter() + .copied() + .fold(OnlineMoments::<3>::default(), OnlineMoments::push); + assert_eq!(merged.count, whole.count); + prec::assert_relative_eq!( + merged.mean().unwrap(), + whole.mean().unwrap(), + max_relative = 1e-14 + ); + prec::assert_relative_eq!( + merged.variance().unwrap(), + whole.variance().unwrap(), + max_relative = 1e-13 + ); + prec::assert_relative_eq!( + merged.skewness().unwrap(), + whole.skewness().unwrap(), + max_relative = 1e-12 + ); + } + } + + #[test] + fn merge_with_empty_is_identity() { + let a = [1.0_f64, 2.0, 3.0] + .iter() + .copied() + .fold(OnlineMoments::<2>::default(), OnlineMoments::push); + let empty = OnlineMoments::<2>::default(); + assert_eq!(a.merge(empty).variance(), Some(1.0)); + let a = [1.0_f64, 2.0, 3.0] + .iter() + .copied() + .fold(OnlineMoments::<2>::default(), OnlineMoments::push); + let empty = OnlineMoments::<2>::default(); + assert_eq!(empty.merge(a).variance(), Some(1.0)); + } + #[test] fn order_3_mean_and_variance_match_order_2() { let data = [2.0_f64, 4.0, 4.0, 4.0, 5.0, 5.0, 7.0, 9.0]; From 5e417980a5216be562346e0e677825b5b8f20b11 Mon Sep 17 00:00:00 2001 From: Felix Agene Date: Sun, 26 Jul 2026 18:44:09 -0500 Subject: [PATCH 2/2] test: build the OnlineMoments offset tests under no_std The tests collected their inputs into `Vec`, which does not exist under `--no-default-features`: the crate is `no_std` without `alloc`. Folding straight off the iterators gives the same accumulations with no allocation; `merge_reconciles_different_offsets` needs the data three times, so it takes closures returning fresh iterators rather than one collected copy. CI did not catch this because the test job runs only default features, and the feature-powerset job passes `--no-dev-deps`, so test code is never compiled without std. Co-Authored-By: Claude Opus 5 (1M context) --- src/statistics/online.rs | 30 +++++++++++------------------- 1 file changed, 11 insertions(+), 19 deletions(-) diff --git a/src/statistics/online.rs b/src/statistics/online.rs index eab11bc1..09938237 100644 --- a/src/statistics/online.rs +++ b/src/statistics/online.rs @@ -300,12 +300,10 @@ mod tests { for exp in [0i32, 10, 20, 30, 40] { let offset = f64::powi(2.0, exp); - let data: Vec = (0..PERIOD * REPEATS) + // Folded straight off the iterator rather than collected, so the + // test also builds under `no_std`, where there is no `Vec`. + let s = (0..PERIOD * REPEATS) .map(|i| offset + (i % PERIOD) as f64 * STEP) - .collect(); - let s = data - .iter() - .copied() .fold(OnlineMoments::<2>::default(), OnlineMoments::push); prec::assert_relative_eq!( s.variance().unwrap(), @@ -328,20 +326,14 @@ mod tests { /// offsets; check it still matches a single chain on offset data. #[test] fn merge_reconciles_different_offsets() { - let a: Vec = (0..500).map(|i| 1e12 + i as f64 * 1e-3).collect(); - let b: Vec = (0..500).map(|i| 1e12 + 5.0 + i as f64 * 1e-3).collect(); - let ma = a - .iter() - .copied() - .fold(OnlineMoments::<2>::default(), OnlineMoments::push); - let mb = b - .iter() - .copied() - .fold(OnlineMoments::<2>::default(), OnlineMoments::push); - let whole: Vec = a.iter().chain(b.iter()).copied().collect(); - let mw = whole - .iter() - .copied() + // Closures returning fresh iterators, so the same data can be folded + // three ways without collecting it (there is no `Vec` under `no_std`). + let a = || (0..500).map(|i| 1e12 + i as f64 * 1e-3); + let b = || (0..500).map(|i| 1e12 + 5.0 + i as f64 * 1e-3); + let ma = a().fold(OnlineMoments::<2>::default(), OnlineMoments::push); + let mb = b().fold(OnlineMoments::<2>::default(), OnlineMoments::push); + let mw = a() + .chain(b()) .fold(OnlineMoments::<2>::default(), OnlineMoments::push); prec::assert_relative_eq!( ma.merge(mb).variance().unwrap(),