hpr_analysis/
statistics.rs

1//! Summaries of a sample of numbers: the spread of a Monte Carlo run's apogees, say.
2//!
3//! A [`Distribution`] keeps every value, sorted, and the number of samples that were tried, so a
4//! sample that failed or never gave a value (a flight with no apogee) is counted, not dropped:
5//! its share is [`Distribution::missing`] of [`Distribution::attempted`], and a probability is
6//! reported as the bounds those unknowns allow ([`Distribution::share_at_least`]).
7//!
8//! - **Mean and standard deviation** are taken on the values shifted by the smallest one, the
9//!   standard deviation by the two-pass formula with `n − 1` (T. F. Chan, G. H. Golub and R. J.
10//!   LeVeque, "Algorithms for computing the sample variance: analysis and recommendations", *The
11//!   American Statistician* 37(3), 242–247, 1983, <https://doi.org/10.2307/2683386>). Shifting
12//!   by a value of the sample keeps the sums small, and makes a sample of equal values give that
13//!   value and a deviation of exactly zero.
14//! - **Quantiles** are Hyndman and Fan's definition 7, linear between order statistics, the
15//!   default of R and NumPy: with the values sorted `x₀ ≤ … ≤ xₙ₋₁` and `h = (n − 1) p`,
16//!   `Q(p) = x⌊h⌋ + (h − ⌊h⌋)(x⌊h⌋₊₁ − x⌊h⌋)` (R. J. Hyndman and Y. Fan, "Sample quantiles in
17//!   statistical packages", *The American Statistician* 50(4), 361–365, 1996,
18//!   <https://doi.org/10.2307/2684934>).
19//!
20//! Every sum runs over the sorted values in order, so a summary is bit-for-bit the same however
21//! the values were computed, in parallel or not.
22
23use serde::{Deserialize, Serialize};
24
25use crate::error::AnalysisError;
26
27/// The values a sample of runs gave, sorted, and how many runs were tried. It serializes as
28/// those two, and reads back through [`Distribution::new`]'s checks.
29#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
30#[serde(try_from = "DistributionData")]
31pub struct Distribution {
32    attempted: usize,
33    sorted: Vec<f64>,
34}
35
36/// The serialized form of a [`Distribution`].
37#[derive(Deserialize)]
38#[serde(deny_unknown_fields)]
39struct DistributionData {
40    attempted: usize,
41    sorted: Vec<f64>,
42}
43
44impl TryFrom<DistributionData> for Distribution {
45    type Error = AnalysisError;
46
47    fn try_from(data: DistributionData) -> Result<Self, AnalysisError> {
48        Self::new(data.sorted, data.attempted)
49    }
50}
51
52/// The usual numbers of a [`Distribution`], for a report.
53#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
54pub struct Summary {
55    /// The runs tried.
56    pub attempted: usize,
57    /// The runs that gave a value.
58    pub count: usize,
59    /// The mean, `None` with no values.
60    pub mean: Option<f64>,
61    /// The sample standard deviation (`n − 1`), `None` with fewer than two values.
62    pub standard_deviation: Option<f64>,
63    /// The smallest value.
64    pub min: Option<f64>,
65    /// The 5th percentile ([`Distribution::quantile`]).
66    pub p05: Option<f64>,
67    /// The median.
68    pub p50: Option<f64>,
69    /// The 95th percentile.
70    pub p95: Option<f64>,
71    /// The largest value.
72    pub max: Option<f64>,
73}
74
75/// Bounds on the share of all the runs tried whose value passed a test: `low` counts a run with
76/// no value as failing it, `high` as passing.
77#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
78pub struct Share {
79    /// The runs with a value that passed, over every run tried.
80    pub low: f64,
81    /// The same, with every run that gave no value counted as passing.
82    pub high: f64,
83}
84
85impl Distribution {
86    /// The distribution of `values`, from `attempted` runs (those that gave no value make up the
87    /// difference).
88    ///
89    /// # Errors
90    ///
91    /// - [`AnalysisError::Domain`] for a value that isn't finite.
92    /// - [`AnalysisError::Count`] for more values than runs.
93    pub fn new(values: Vec<f64>, attempted: usize) -> Result<Self, AnalysisError> {
94        if values.len() > attempted {
95            return Err(AnalysisError::Count {
96                what: "values in a distribution, against the runs tried",
97                count: values.len(),
98                limit: attempted,
99            });
100        }
101        if let Some(&bad) = values.iter().find(|v| !v.is_finite()) {
102            return Err(AnalysisError::Domain {
103                what: "value in a distribution",
104                value: bad,
105            });
106        }
107        let mut sorted = values;
108        sorted.sort_by(f64::total_cmp);
109        Ok(Self { attempted, sorted })
110    }
111
112    /// The runs tried.
113    pub fn attempted(&self) -> usize {
114        self.attempted
115    }
116
117    /// The runs that gave a value.
118    pub fn count(&self) -> usize {
119        self.sorted.len()
120    }
121
122    /// The runs that gave no value: failed, or without the quantity asked for.
123    pub fn missing(&self) -> usize {
124        self.attempted - self.sorted.len()
125    }
126
127    /// The values, smallest first.
128    pub fn sorted(&self) -> &[f64] {
129        &self.sorted
130    }
131
132    /// The smallest value, `None` with no values.
133    pub fn min(&self) -> Option<f64> {
134        self.sorted.first().copied()
135    }
136
137    /// The largest value, `None` with no values.
138    pub fn max(&self) -> Option<f64> {
139        self.sorted.last().copied()
140    }
141
142    /// The mean, `x₀ + Σ(xᵢ − x₀)/n` with `x₀` the smallest value; `None` with no values.
143    pub fn mean(&self) -> Option<f64> {
144        let shift = self.min()?;
145        Some(shift + self.shifted_mean(shift))
146    }
147
148    /// The sample standard deviation, `√(Σ(dᵢ − d̄)²/(n − 1))` with `dᵢ = xᵢ − x₀`; `None` with
149    /// fewer than two values.
150    pub fn standard_deviation(&self) -> Option<f64> {
151        if self.sorted.len() < 2 {
152            return None;
153        }
154        let shift = self.min()?;
155        let mean = self.shifted_mean(shift);
156        let sum_squares: f64 = self
157            .sorted
158            .iter()
159            .map(|&x| {
160                let d = (x - shift) - mean;
161                d * d
162            })
163            .sum();
164        // Cast: a count of values is far below 2⁵³.
165        Some((sum_squares / (self.sorted.len() - 1) as f64).sqrt())
166    }
167
168    /// The mean of the values less `shift`.
169    fn shifted_mean(&self, shift: f64) -> f64 {
170        let sum: f64 = self.sorted.iter().map(|&x| x - shift).sum();
171        // Cast: a count of values is far below 2⁵³.
172        sum / self.sorted.len() as f64
173    }
174
175    /// The `p` quantile by Hyndman and Fan's definition 7 (the module's docs); `None` with no
176    /// values.
177    ///
178    /// # Errors
179    ///
180    /// [`AnalysisError::Domain`] for `p` outside `[0, 1]`.
181    pub fn quantile(&self, p: f64) -> Result<Option<f64>, AnalysisError> {
182        if !(0.0..=1.0).contains(&p) {
183            return Err(AnalysisError::Domain {
184                what: "quantile probability",
185                value: p,
186            });
187        }
188        let Some(&last) = self.sorted.last() else {
189            return Ok(None);
190        };
191        // Cast: a count of values is far below 2⁵³, and `h` lies in `[0, n − 1]`.
192        let h = (self.sorted.len() - 1) as f64 * p;
193        let below = h.floor();
194        let index = below as usize;
195        let Some(&above) = self.sorted.get(index + 1) else {
196            return Ok(Some(last));
197        };
198        let low = self.sorted[index];
199        Ok(Some(low + (h - below) * (above - low)))
200    }
201
202    /// Bounds on the share of the runs tried whose value is at least `threshold` ([`Share`]);
203    /// `None` with no runs tried, when there is no share to bound. An infinite threshold is
204    /// allowed.
205    ///
206    /// # Errors
207    ///
208    /// [`AnalysisError::Domain`] for a threshold that is not a number.
209    pub fn share_at_least(&self, threshold: f64) -> Result<Option<Share>, AnalysisError> {
210        if threshold.is_nan() {
211            return Err(AnalysisError::Domain {
212                what: "threshold of a share",
213                value: threshold,
214            });
215        }
216        if self.attempted == 0 {
217            return Ok(None);
218        }
219        let passed = self.sorted.len() - self.sorted.partition_point(|&x| x < threshold);
220        // Cast: counts far below 2⁵³.
221        let attempted = self.attempted as f64;
222        Ok(Some(Share {
223            low: passed as f64 / attempted,
224            high: (passed + self.missing()) as f64 / attempted,
225        }))
226    }
227
228    /// The usual numbers, for a report: [`Summary`].
229    pub fn summary(&self) -> Summary {
230        // The probabilities are in [0, 1].
231        let quantile = |p: f64| self.quantile(p).ok().flatten();
232        Summary {
233            attempted: self.attempted,
234            count: self.count(),
235            mean: self.mean(),
236            standard_deviation: self.standard_deviation(),
237            min: self.min(),
238            p05: quantile(0.05),
239            p50: quantile(0.5),
240            p95: quantile(0.95),
241            max: self.max(),
242        }
243    }
244}
245
246#[cfg(test)]
247mod tests {
248    use super::*;
249
250    #[test]
251    fn equal_values_give_their_value_and_no_spread() {
252        // A sum of 0.1s divided by their count isn't 0.1 in floating point; shifted by the
253        // smallest value, it is.
254        let d = Distribution::new(vec![0.1; 7], 7).unwrap();
255        assert_eq!(d.mean(), Some(0.1));
256        assert_eq!(d.standard_deviation(), Some(0.0));
257        assert_eq!(d.quantile(0.37).unwrap(), Some(0.1));
258    }
259
260    #[test]
261    fn moments_and_quantiles_by_hand() {
262        let d = Distribution::new(vec![4.0, 1.0, 3.0, 2.0], 6).unwrap();
263        assert_eq!(d.sorted(), &[1.0, 2.0, 3.0, 4.0]);
264        assert_eq!((d.count(), d.missing(), d.attempted()), (4, 2, 6));
265        assert_eq!(d.mean(), Some(2.5));
266        // Σ(x − 2.5)² = 5, over n − 1 = 3.
267        assert_eq!(d.standard_deviation(), Some((5.0_f64 / 3.0).sqrt()));
268        // h = 3p: p = 0.5 is halfway from 2 to 3; p = 0.9 is 0.7 of the way from 3 to 4. NumPy's
269        // `np.quantile([1, 2, 3, 4], [0, 0.5, 0.9, 1])` gives the same: 1, 2.5, 3.7, 4.
270        assert_eq!(d.quantile(0.0).unwrap(), Some(1.0));
271        assert_eq!(d.quantile(0.5).unwrap(), Some(2.5));
272        assert!((d.quantile(0.9).unwrap().unwrap() - 3.7).abs() < 1e-15);
273        assert_eq!(d.quantile(1.0).unwrap(), Some(4.0));
274        let summary = d.summary();
275        assert_eq!((summary.min, summary.max), (Some(1.0), Some(4.0)));
276        assert_eq!(summary.p50, Some(2.5));
277    }
278
279    #[test]
280    fn a_share_is_bounded_by_the_runs_with_no_value() {
281        // Two of six runs gave nothing: at least 2/6 and at most 4/6 reached 3.
282        let d = Distribution::new(vec![1.0, 2.0, 3.0, 4.0], 6).unwrap();
283        let share = d.share_at_least(3.0).unwrap().unwrap();
284        assert_eq!((share.low, share.high), (2.0 / 6.0, 4.0 / 6.0));
285        let all = Distribution::new(vec![1.0, 2.0], 2).unwrap();
286        let share = all.share_at_least(1.0).unwrap().unwrap();
287        assert_eq!((share.low, share.high), (1.0, 1.0));
288        let share = all.share_at_least(f64::INFINITY).unwrap().unwrap();
289        assert_eq!((share.low, share.high), (0.0, 0.0));
290        let share = all.share_at_least(f64::NEG_INFINITY).unwrap().unwrap();
291        assert_eq!((share.low, share.high), (1.0, 1.0));
292        // No runs, no share; a threshold that isn't a number is refused.
293        let none = Distribution::new(vec![], 0).unwrap();
294        assert_eq!(none.share_at_least(1.0).unwrap(), None);
295        assert!(matches!(
296            all.share_at_least(f64::NAN),
297            Err(AnalysisError::Domain {
298                what: "threshold of a share",
299                ..
300            })
301        ));
302    }
303
304    #[test]
305    fn empty_and_single_samples() {
306        let empty = Distribution::new(vec![], 3).unwrap();
307        assert_eq!(empty.mean(), None);
308        assert_eq!(empty.quantile(0.5).unwrap(), None);
309        assert_eq!(empty.missing(), 3);
310        let one = Distribution::new(vec![2.0], 1).unwrap();
311        assert_eq!(one.mean(), Some(2.0));
312        assert_eq!(one.standard_deviation(), None);
313        assert_eq!(one.quantile(0.3).unwrap(), Some(2.0));
314    }
315
316    #[test]
317    fn a_distribution_reads_back_through_its_checks() {
318        let d = Distribution::new(vec![3.0, 1.0, 2.0], 4).unwrap();
319        let json = serde_json::to_string(&d).unwrap();
320        assert_eq!(json, r#"{"attempted":4,"sorted":[1.0,2.0,3.0]}"#);
321        assert_eq!(serde_json::from_str::<Distribution>(&json).unwrap(), d);
322        // Unsorted values are sorted on the way in; more values than runs are refused.
323        let unsorted: Distribution =
324            serde_json::from_str(r#"{"attempted":3,"sorted":[3.0,1.0,2.0]}"#).unwrap();
325        assert_eq!(unsorted.sorted(), &[1.0, 2.0, 3.0]);
326        let error = serde_json::from_str::<Distribution>(r#"{"attempted":1,"sorted":[3.0,1.0]}"#)
327            .unwrap_err();
328        assert!(error.to_string().contains("more than 1"), "{error}");
329    }
330
331    #[test]
332    fn bad_inputs_are_refused() {
333        assert!(matches!(
334            Distribution::new(vec![1.0, f64::NAN], 2),
335            Err(AnalysisError::Domain { what: "value in a distribution", value }) if value.is_nan()
336        ));
337        assert!(matches!(
338            Distribution::new(vec![1.0, 2.0], 1),
339            Err(AnalysisError::Count {
340                count: 2,
341                limit: 1,
342                ..
343            })
344        ));
345        let d = Distribution::new(vec![1.0], 1).unwrap();
346        for p in [-0.1, 1.1, f64::NAN] {
347            assert!(matches!(
348                d.quantile(p),
349                Err(AnalysisError::Domain { what: "quantile probability", value })
350                    if value.to_bits() == p.to_bits()
351            ));
352        }
353    }
354}
hpr_analysis/statistics.rs

hpr_analysis/
statistics.rs
