shape_runtime/simd_rolling/
std.rs

1//! Rolling standard deviation using Welford's algorithm
2
3/// Rolling std (auto-selects Welford algorithm, always optimal)
4#[inline]
5pub fn rolling_std(data: &[f64], window: usize) -> Vec<f64> {
6    rolling_std_welford(data, window) // Welford is always best
7}
8
9/// Welford's algorithm for numerically stable rolling standard deviation
10///
11/// This is significantly faster than the naive two-pass approach and
12/// more numerically stable than the sum-of-squares method.
13///
14/// # Performance
15/// - Expected speedup: 10-50x over naive O(n*w) two-pass
16/// - Complexity: O(n)
17/// - Numerical stability: Superior to sum-of-squares
18pub fn rolling_std_welford(data: &[f64], window: usize) -> Vec<f64> {
19    let n = data.len();
20
21    if window == 0 || window > n {
22        return vec![f64::NAN; n];
23    }
24
25    let mut result = vec![f64::NAN; n];
26
27    // Initialize with first window using Welford's algorithm
28    let mut mean = 0.0;
29    let mut m2 = 0.0; // Sum of squared differences from mean
30
31    for i in 0..window {
32        let delta = data[i] - mean;
33        mean += delta / (i + 1) as f64;
34        let delta2 = data[i] - mean;
35        m2 += delta * delta2;
36    }
37
38    result[window - 1] = (m2 / window as f64).sqrt();
39
40    // Slide window: remove oldest, add newest
41    for i in window..n {
42        let old_val = data[i - window];
43        let new_val = data[i];
44
45        // Remove old value from running statistics
46        let old_mean = mean;
47        mean -= (old_val - mean) / window as f64;
48        m2 -= (old_val - mean) * (old_val - old_mean);
49
50        // Add new value
51        let delta = new_val - mean;
52        mean += delta / window as f64;
53        let delta2 = new_val - mean;
54        m2 += delta * delta2;
55
56        // Ensure non-negative due to floating point errors
57        result[i] = (m2.max(0.0) / window as f64).sqrt();
58    }
59
60    result
61}
62
63#[cfg(test)]
64mod tests {
65    use super::*;
66
67    #[test]
68    fn test_rolling_std_welford() {
69        let data = vec![2.0, 4.0, 4.0, 4.0, 5.0, 5.0, 7.0, 9.0];
70        let result = rolling_std_welford(&data, 4);
71
72        // First 3 are NaN
73        for i in 0..3 {
74            assert!(result[i].is_nan());
75        }
76
77        // Check that we get reasonable std values
78        assert!(result[3] > 0.0 && result[3] < 2.0); // std of [2,4,4,4]
79        assert!(result[7] > 0.0 && result[7] < 3.0); // std of [5,5,7,9]
80    }
81}
shape_runtime/simd_rolling/std.rs

shape_runtime/simd_rolling/
std.rs
