online moments

[YeungOnion]

• feat: add OnlineMoments accumulator using Welford
• fix: wrap online stats with iter_statistics for Numerical Stability of Variance #376
• fix: drop excess precision requirements

[codecov]

Codecov Report

❌ Patch coverage is 80.10204% with 39 lines in your changes missing coverage. Please review.
✅ Project coverage is 94.67%. Comparing base (14d3f8f) to head (faaf80a).
⚠️ Report is 1 commits behind head on main.

Files with missing lines	Patch %	Lines
src/statistics/accumulate.rs	11.53%	23 Missing ⚠️
src/statistics/online.rs	89.47%	16 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #394      +/-   ##
==========================================
- Coverage   94.69%   94.67%   -0.02%     
==========================================
  Files          59       61       +2     
  Lines       13068    13322     +254     
==========================================
+ Hits        12375    12613     +238     
- Misses        693      709      +16     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[YeungOnion]

forgot to mention that this was likely only on my mind due to #324, so thanks for the feedback there.

[Evian-Zhang]

Wow wonderful work! I can't wait to use this feature!

However, I have two questions:

1. Why OnlineMean is OnlineMoments<2> instead of OnlineMoments<1>? I'm not expert in statistics, is this a typo or intentional use?

2. It seems to me that push is actually called several times for each element since

   statrs/src/statistics/accumulate.rs

   Lines 53 to 61 in 10cf6d6

   	fn push(self, x: f64) -> Self {
   	(
   	self.0.push(x),
   	self.1.push(x),
   	self.2.push(x),
   	self.3.push(x),
   	self.4.push(x),
   	)
   	}

   . Moreover, if we write

   let (mean, var, skewness): (......) = data.iter().copied().fold(Default::default(), Accumulate::push);

   Then the mean value can be calculated by mean.mean(), var.mean() and skewness.mean(), which is a little bit redundant to me.

Based on your code, it occurs to me that we can really write declarative code like this:

    let (OnlineMean(mean)) = data.iter().copied().fold(Default::default(), Accumulate::push).get();
    // or
    let (OnlineMean(mean), OnlineVariance(variance)) = data.iter().copied().fold(Default::default(), Accumulate::push).get();
    // or
    let (OnlineVariance(variance), OnlineSkewness(skewness)) = data.iter().copied().fold(Default::default(), Accumulate::push).get();
    // then mean, variance, and skewness is Option<f64> directly

if we use the following type magic:

pub trait Moment {
    const ORDER: usize;
    // This method is not meant to be called from user, so we can
    // have less type restrictions. In theory, we can even restrict
    // MS only to the tuple types that include Self.
    fn from<MS: Moments>(acc: &Accumulate<MS>) -> Self;
}
pub trait Moments: Sized {
    // This will always be compiled into constants.
    // In the future, this can be directly an associated const.
    fn order() -> usize;
    fn from(acc: &Accumulate<Self>) -> Self;
}

impl<M1: Moment, M2: Moment> Moments for (M1, M2) {
    fn order() -> usize {
        std::cmp::max(M1::ORDER, M2::ORDER)
    }

    fn from(acc: &Accumulate<Self>) -> Self {
        (M1::from(acc), M2::from(acc))
    }
}
// Also implements for (M1), (M1, M2, M3), ...

pub struct Accumulate<MS: Moments> {
    count: u64,
    // Some type magics can be done here to make
    // the array length equal to MS::order(). But since
    // we only need up to three-order moment, a simple
    // [_; 3] is enough.
    m: [f64; 3],
    phantom: core::marker::PhantomData<MS>,
}

impl<MS: Moments> Default for Accumulate<MS> {
    fn default() -> Self {
        Self {
            count: 0,
            m: [0.0; 3],
            phantom: core::marker::PhantomData::default()
        }
    }
}

impl<MS: Moments> Accumulate<MS> {
    pub fn push(mut self, x: f64) -> Self {
        self.count += 1;
        let n = self.count as f64;

        // Welford / Pebay (2008) central moment update. Update order: M3
        // before M2 before mean; each step uses the previous observation's
        // lower-order accumulators.
        let delta = x - self.m[0];
        let delta_n = delta / n;
        let new_mean = self.m[0] + delta_n;
        let delta2 = x - new_mean;

        if MS::order() >=2  {
            let old_m2 = self.m[1];
            if MS::order() >=3 {
                let inc = delta * (delta_n * delta_n) * (n - 1.0) * (n - 2.0) - 3.0 * delta_n * old_m2;
                self.m[2] += inc;
            }
            self.m[1] += delta * delta2;
        }

        self.m[0] = new_mean;
        self
    }

    pub fn get(&self) -> MS {
        MS::from(self)
    }
}

pub struct OnlineMean(pub Option<f64>);
impl Moment for OnlineMean {
    const ORDER: usize = 1;
    fn from<MS: Moments>(acc: &Accumulate<MS>) -> Self {
        if acc.count == 0 {
            Self(None)
        } else {
            Self(Some(acc.m[0]))
        }
    }
}

pub struct OnlineVariance(pub Option<f64>);
impl Moment for OnlineVariance {
    const ORDER: usize = 2;
    fn from<MS: Moments>(acc: &Accumulate<MS>) -> Self {
        if acc.count < 2 {
            Self(None)
        } else {
            Self(Some(acc.m[1] / (acc.count - 1) as f64))
        }
    }
}

// Also implements for OlineSkewness, OnlineStdDev, OnlinePopulationVariance...

The code snippet above has been tested by me that is consistent to your implementation. The difference is that the usage code can be more declarative, and push is only called one time for each element.

[Evian-Zhang]

@YeungOnion What do you think of this idea? I'm happy to send a PR after #408 is merged (otherwise that PR has to be refactored a lot).

[YeungOnion]

yeah, something like that could be nice, I do think there's a way to surface that welford provides all lower order moments and doesn't need separate reducer states. I just merged #408 so you can take a stab at it. I do think part of it was my poor choice in names of type aliases for OnlineMoments. I do like the destructure pattern instead of duplicating a call to mean method. But the .get finalize seems like a great call for this approach instead of delegating finalization to .mean