rustyqlib/core/montecarlo/
stats.rs

1//! Simulation statistics: mean / standard error from accumulated sums
2//! (the shape parallel path loops naturally produce) and a numerically
3//! stable Welford accumulator for streaming use.
4
5/// A simulation estimate: sample mean, its standard error, and the
6/// sample size.
7#[derive(Debug, Clone, Copy)]
8pub struct SimStats {
9    pub mean: f64,
10    pub std_err: f64,
11    pub n: usize,
12}
13
14/// Mean and standard error of the mean from `sum` and `sum_sq` of `n`
15/// samples — the reduction shape of a parallel path loop.
16pub fn mean_std_err(sum: f64, sum_sq: f64, n: usize) -> (f64, f64) {
17    let nf = n as f64;
18    let mean = sum / nf;
19    let var = (sum_sq / nf - mean * mean).max(0.0);
20    (mean, (var / nf).sqrt())
21}
22
23/// Welford's online mean/variance accumulator: numerically stable for
24/// long streams and mergeable across partial results.
25#[derive(Debug, Clone, Copy, Default)]
26pub struct RunningStats {
27    n: usize,
28    mean: f64,
29    m2: f64,
30}
31
32impl RunningStats {
33    pub fn new() -> Self {
34        Self::default()
35    }
36
37    pub fn push(&mut self, x: f64) {
38        self.n += 1;
39        let delta = x - self.mean;
40        self.mean += delta / self.n as f64;
41        self.m2 += delta * (x - self.mean);
42    }
43
44    /// Merge another accumulator (parallel reduction).
45    pub fn merge(&mut self, other: &RunningStats) {
46        if other.n == 0 {
47            return;
48        }
49        if self.n == 0 {
50            *self = *other;
51            return;
52        }
53        let n1 = self.n as f64;
54        let n2 = other.n as f64;
55        let delta = other.mean - self.mean;
56        let n = n1 + n2;
57        self.mean += delta * n2 / n;
58        self.m2 += other.m2 + delta * delta * n1 * n2 / n;
59        self.n += other.n;
60    }
61
62    pub fn count(&self) -> usize {
63        self.n
64    }
65
66    pub fn mean(&self) -> f64 {
67        self.mean
68    }
69
70    /// Population variance (divide by n, matching [`mean_std_err`]).
71    pub fn variance(&self) -> f64 {
72        if self.n == 0 { 0.0 } else { (self.m2 / self.n as f64).max(0.0) }
73    }
74
75    pub fn std_err(&self) -> f64 {
76        if self.n == 0 { 0.0 } else { (self.variance() / self.n as f64).sqrt() }
77    }
78
79    pub fn stats(&self) -> SimStats {
80        SimStats { mean: self.mean(), std_err: self.std_err(), n: self.n }
81    }
82}
83
84#[cfg(test)]
85mod tests {
86    use super::*;
87
88    #[test]
89    fn welford_matches_the_two_pass_computation() {
90        let xs: Vec<f64> = (0..1000).map(|i| ((i * 37) % 101) as f64 * 0.13 - 5.0).collect();
91        let mut running = RunningStats::new();
92        for &x in &xs {
93            running.push(x);
94        }
95        let sum: f64 = xs.iter().sum();
96        let sum_sq: f64 = xs.iter().map(|x| x * x).sum();
97        let (mean, se) = mean_std_err(sum, sum_sq, xs.len());
98        assert!((running.mean() - mean).abs() < 1e-10);
99        assert!((running.std_err() - se).abs() < 1e-10);
100    }
101
102    #[test]
103    fn merged_accumulators_equal_a_single_pass() {
104        let xs: Vec<f64> = (0..500).map(|i| (i as f64).sin() * 3.0).collect();
105        let mut whole = RunningStats::new();
106        let mut left = RunningStats::new();
107        let mut right = RunningStats::new();
108        for (i, &x) in xs.iter().enumerate() {
109            whole.push(x);
110            if i < 200 { left.push(x) } else { right.push(x) }
111        }
112        left.merge(&right);
113        assert_eq!(left.count(), whole.count());
114        assert!((left.mean() - whole.mean()).abs() < 1e-12);
115        assert!((left.variance() - whole.variance()).abs() < 1e-12);
116    }
117}
rustyqlib/core/montecarlo/stats.rs

rustyqlib/core/montecarlo/
stats.rs
