lindhard/tally/moments.rs
1//! One-pass, mergeable central moments with standard errors.
2//!
3//! # Accumulator
4//!
5//! [`Moments`] keeps the count `n`, the mean, and the centred power sums
6//! `M_p = Σ (x_i - mean)^p` for `p = 2..=8`. Two accumulators `A`, `B` combine
7//! with the arbitrary-order pairwise update of P. Pébay, "Formulas for robust,
8//! one-pass parallel computation of covariances and arbitrary-order
9//! statistical moments", Sandia report SAND2008-6212 (2008), which generalises
10//! the second-order update of T. F. Chan, G. H. Golub and R. J. LeVeque,
11//! "Updating formulae and a pairwise algorithm for computing sample
12//! variances", Stanford report STAN-CS-79-773 (1979). With `δ = mean_B -
13//! mean_A` and `n = n_A + n_B`,
14//!
15//! ```text
16//! M_p = M_p^A + M_p^B
17//! + Σ_{k=1}^{p-2} C(p,k) δ^k [ (-n_B/n)^k M_{p-k}^A + (n_A/n)^k M_{p-k}^B ]
18//! + (n_A n_B δ / n)^p [ n_B^{1-p} - (-n_A)^{1-p} ].
19//! ```
20//!
21//! Adding one sample is a merge with `n_B = 1`, `M^B = 0`. The update is exact
22//! in real arithmetic; in floating point the result depends on the order of
23//! merges, so it is **deterministic for a fixed order** (bit-identical) and
24//! agrees across different splittings to rounding. [`crate::rng::run_particles`]
25//! merges in chunk order, which fixes the order independently of the thread
26//! count.
27//!
28//! # Conventions
29//!
30//! Central moments are the method-of-moments estimates `μ_p = M_p / n` (no
31//! `n - 1` correction), as used for fitting Pearson distributions:
32//! standard deviation `σ = μ_2^(1/2)`, skewness `γ = μ_3 / μ_2^(3/2)`,
33//! kurtosis `β = μ_4 / μ_2^2` (**not** excess kurtosis: a Gaussian has
34//! `β = 3`). For a depth distribution these are `Rp`, `ΔRp`, `γ`, `β`.
35//!
36//! # Standard errors
37//!
38//! Large-sample standard errors from the asymptotic covariance of sample
39//! central moments, H. Cramér, *Mathematical Methods of Statistics*
40//! (Princeton, 1946), §27.7:
41//!
42//! ```text
43//! n Cov(m_j, m_k) = μ_{j+k} - μ_j μ_k - j μ_{j-1} μ_{k+1} - k μ_{j+1} μ_{k-1}
44//! + j k μ_{j-1} μ_{k-1} μ_2,
45//! ```
46//!
47//! propagated to `σ`, `γ` and `β` by the delta method (first-order Taylor
48//! expansion), with the population moments replaced by their estimates up to
49//! `μ_8`. The standard error of the mean is `σ / n^(1/2)`. No normality is
50//! assumed; for a Gaussian these reduce to the familiar `(6/n)^(1/2)` for `γ`
51//! and `(24/n)^(1/2)` for `β` (checked in a test). They are asymptotic, so
52//! treat them as approximate for small `n` or heavy tails.
53
54use serde::{Deserialize, Serialize};
55
56/// Highest power sum kept.
57const MAX_ORDER: usize = 8;
58
59/// Mergeable accumulator of the count, mean and centred power sums up to
60/// order 8. See the module docs for the update formula.
61#[derive(Debug, Clone, PartialEq, Default, Serialize, Deserialize)]
62pub struct Moments {
63 /// Number of samples.
64 pub n: u64,
65 /// Running mean.
66 pub mean: f64,
67 /// Centred power sums `M_2 .. M_8` (index `p - 2`).
68 pub power_sums: [f64; MAX_ORDER - 1],
69}
70
71fn binom(n: usize, k: usize) -> f64 {
72 (0..k).fold(1.0, |acc, i| acc * (n - i) as f64 / (i + 1) as f64)
73}
74
75impl Moments {
76 /// Empty accumulator.
77 pub fn new() -> Self {
78 Self::default()
79 }
80
81 /// `M_p` for `p` in `0..=8` (`M_0 = n`, `M_1 = 0`).
82 fn m(&self, p: usize) -> f64 {
83 match p {
84 0 => self.n as f64,
85 1 => 0.0,
86 _ => self.power_sums[p - 2],
87 }
88 }
89
90 /// Add one sample.
91 pub fn push(&mut self, x: f64) {
92 self.merge_parts(1, x, &[0.0; MAX_ORDER - 1]);
93 }
94
95 /// Fold `other` into `self` (Pébay 2008).
96 pub fn merge(&mut self, other: &Moments) {
97 self.merge_parts(other.n, other.mean, &other.power_sums);
98 }
99
100 fn merge_parts(&mut self, n_b: u64, mean_b: f64, sums_b: &[f64; MAX_ORDER - 1]) {
101 if n_b == 0 {
102 return;
103 }
104 if self.n == 0 {
105 self.n = n_b;
106 self.mean = mean_b;
107 self.power_sums = *sums_b;
108 return;
109 }
110 let (na, nb) = (self.n as f64, n_b as f64);
111 let n = na + nb;
112 let delta = mean_b - self.mean;
113 let m_b = |p: usize| match p {
114 0 => nb,
115 1 => 0.0,
116 _ => sums_b[p - 2],
117 };
118 let mut out = [0.0; MAX_ORDER - 1];
119 for (idx, slot) in out.iter_mut().enumerate() {
120 let p = idx + 2;
121 let mut s = self.m(p) + m_b(p);
122 for k in 1..=p - 2 {
123 s += binom(p, k)
124 * delta.powi(k as i32)
125 * ((-nb / n).powi(k as i32) * self.m(p - k)
126 + (na / n).powi(k as i32) * m_b(p - k));
127 }
128 let pi = p as i32;
129 s += (na * nb * delta / n).powi(pi) * (nb.powi(1 - pi) - (-na).powi(1 - pi));
130 *slot = s;
131 }
132 self.power_sums = out;
133 self.mean += delta * nb / n;
134 self.n += n_b;
135 }
136
137 /// Central moment estimate `μ_p = M_p / n` (`μ_0 = 1`, `μ_1 = 0`).
138 pub fn central(&self, p: usize) -> f64 {
139 assert!(p <= MAX_ORDER, "central moments are kept up to order 8");
140 match p {
141 0 => 1.0,
142 1 => 0.0,
143 _ => self.power_sums[p - 2] / self.n as f64,
144 }
145 }
146
147 /// Mean, standard deviation, skewness, kurtosis and their standard
148 /// errors. `None` if fewer than two samples or zero variance.
149 pub fn summary(&self) -> Option<MomentSummary> {
150 if self.n < 2 {
151 return None;
152 }
153 let mu = |p: usize| self.central(p);
154 let (m2, m3, m4) = (mu(2), mu(3), mu(4));
155 if m2.is_nan() || m2 <= 0.0 {
156 return None;
157 }
158 let n = self.n as f64;
159 // n Cov(m_j, m_k), Cramér (1946) §27.7.
160 let cov = |j: usize, k: usize| {
161 let (jf, kf) = (j as f64, k as f64);
162 mu(j + k) - mu(j) * mu(k) - jf * mu(j - 1) * mu(k + 1) - kf * mu(j + 1) * mu(k - 1)
163 + jf * kf * mu(j - 1) * mu(k - 1) * m2
164 };
165 let quad = |g: [(usize, f64); 2]| {
166 let mut v = 0.0;
167 for &(j, gj) in &g {
168 for &(k, gk) in &g {
169 v += gj * gk * cov(j, k);
170 }
171 }
172 (v.max(0.0) / n).sqrt()
173 };
174 let sd = m2.sqrt();
175 let skew = m3 / (m2 * sd);
176 let kurt = m4 / (m2 * m2);
177 Some(MomentSummary {
178 n: self.n,
179 mean: self.mean,
180 std_dev: sd,
181 skewness: skew,
182 kurtosis: kurt,
183 mean_std_err: (m2 / n).sqrt(),
184 std_dev_std_err: (cov(2, 2).max(0.0) / n).sqrt() / (2.0 * sd),
185 skewness_std_err: quad([(2, -1.5 * m3 / (m2 * m2 * sd)), (3, 1.0 / (m2 * sd))]),
186 kurtosis_std_err: quad([(2, -2.0 * m4 / (m2 * m2 * m2)), (4, 1.0 / (m2 * m2))]),
187 })
188 }
189}
190
191/// The first four moments of a sample, with standard errors.
192///
193/// Units: `mean`, `std_dev` and their errors are in the unit of the samples
194/// (metres for depths and lateral positions); `skewness` and `kurtosis` are
195/// dimensionless. For a depth distribution `mean` is the projected range
196/// `Rp`, `std_dev` the straggle `ΔRp`, `skewness` `γ` and `kurtosis` `β`
197/// (Gaussian `β = 3`). See the module docs for conventions.
198#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
199pub struct MomentSummary {
200 /// Number of samples.
201 pub n: u64,
202 /// Mean.
203 pub mean: f64,
204 /// Standard deviation `μ_2^(1/2)`.
205 pub std_dev: f64,
206 /// Skewness `γ = μ_3 / μ_2^(3/2)`.
207 pub skewness: f64,
208 /// Kurtosis `β = μ_4 / μ_2^2` (not excess).
209 pub kurtosis: f64,
210 /// Standard error of `mean`.
211 pub mean_std_err: f64,
212 /// Standard error of `std_dev`.
213 pub std_dev_std_err: f64,
214 /// Standard error of `skewness`.
215 pub skewness_std_err: f64,
216 /// Standard error of `kurtosis`.
217 pub kurtosis_std_err: f64,
218}
219
220#[cfg(test)]
221mod tests {
222 use super::*;
223
224 fn two_pass(xs: &[f64]) -> (f64, Vec<f64>) {
225 let n = xs.len() as f64;
226 let mean = xs.iter().sum::<f64>() / n;
227 let sums = (2..=8)
228 .map(|p| xs.iter().map(|x| (x - mean).powi(p)).sum::<f64>())
229 .collect();
230 (mean, sums)
231 }
232
233 #[test]
234 fn push_matches_two_pass() {
235 // A skewed deterministic sample.
236 let xs: Vec<f64> = (0..500)
237 .map(|i| ((i as f64) * 0.37).sin().exp() * 3.0 + 1.0)
238 .collect();
239 let mut m = Moments::new();
240 xs.iter().for_each(|&x| m.push(x));
241 let (mean, sums) = two_pass(&xs);
242 assert!((m.mean - mean).abs() < 1e-13 * mean.abs());
243 for (a, b) in m.power_sums.iter().zip(&sums) {
244 assert!((a - b).abs() <= 1e-10 * b.abs().max(1e-12), "{a} vs {b}");
245 }
246 }
247
248 #[test]
249 fn merge_of_empty_is_identity() {
250 let mut a = Moments::new();
251 a.push(1.0);
252 a.push(4.0);
253 let before = a.clone();
254 a.merge(&Moments::new());
255 assert_eq!(a, before);
256 let mut e = Moments::new();
257 e.merge(&before);
258 assert_eq!(e, before);
259 }
260
261 #[test]
262 fn normal_theory_standard_errors() {
263 // Exact Gaussian moments with sigma = 1: mu4 = 3, mu6 = 15, mu8 = 105.
264 let m = Moments {
265 n: 100,
266 mean: 0.0,
267 power_sums: [100.0, 0.0, 300.0, 0.0, 1500.0, 0.0, 10500.0],
268 };
269 let s = m.summary().unwrap();
270 assert!((s.skewness_std_err - (6.0f64 / 100.0).sqrt()).abs() < 1e-14);
271 assert!((s.kurtosis_std_err - (24.0f64 / 100.0).sqrt()).abs() < 1e-14);
272 assert!((s.std_dev_std_err - (0.5f64 / 100.0).sqrt()).abs() < 1e-14);
273 assert!((s.mean_std_err - 0.1).abs() < 1e-15);
274 }
275
276 #[test]
277 fn degenerate_samples_have_no_summary() {
278 let mut m = Moments::new();
279 assert!(m.summary().is_none());
280 m.push(2.0);
281 assert!(m.summary().is_none());
282 m.push(2.0);
283 assert!(m.summary().is_none());
284 }
285}