Skip to main content

optirs_nas/evaluation/
statistical.rs

1//! Statistical analysis for evaluation results
2//!
3//! Provides statistical tests, analysis methods, and confidence intervals.
4
5use scirs2_core::numeric::Float;
6use std::cmp::Ordering;
7use std::collections::HashMap;
8use std::fmt::Debug;
9
10use super::benchmark::TestResult;
11use super::types::*;
12use crate::error::{OptimError, Result};
13use crate::numeric::{count_as, scalar_as};
14use crate::EvaluationMetric;
15
16/// Statistical analyzer for evaluation results
17///
18/// `analysis_methods` and `multiple_comparison` were configuration the analyzer
19/// stored and never consulted: it always produced descriptive statistics and
20/// never corrected a p-value, whatever they said. Both now drive behaviour —
21/// see [`StatisticalAnalyzer::analyze`] and
22/// [`StatisticalAnalyzer::adjusted_p_values`] — and the recorded
23/// [`StatisticalTest`]s they operate on are readable through
24/// [`StatisticalAnalyzer::tests`].
25#[derive(Debug)]
26pub struct StatisticalAnalyzer<T: Float + Debug + Send + Sync + 'static> {
27    /// Statistical tests
28    statistical_tests: Vec<StatisticalTest<T>>,
29
30    /// Analyses [`StatisticalAnalyzer::analyze`] performs, in order.
31    analysis_methods: Vec<AnalysisMethod>,
32
33    /// Significance thresholds
34    significance_thresholds: HashMap<String, T>,
35
36    /// Family-wise / false-discovery correction applied by
37    /// [`StatisticalAnalyzer::adjusted_p_values`].
38    multiple_comparison: MultipleComparisonCorrection,
39}
40
41/// Result of running the analyzer's enabled analyses over one sample.
42#[derive(Debug, Clone)]
43pub struct AnalysisReport<T: Float + Debug + Send + Sync + 'static> {
44    /// Descriptive statistics, present when
45    /// [`AnalysisMethod::DescriptiveStatistics`] is enabled.
46    pub descriptive: Option<DescriptiveStats<T>>,
47
48    /// Lag-1 autocorrelation of the sample, present when
49    /// [`AnalysisMethod::CorrelationAnalysis`] is enabled and the sample has at
50    /// least two points and non-zero variance. `None` for a constant or
51    /// too-short sample: an autocorrelation is undefined there, and reporting
52    /// `0` would claim "no correlation" where the truth is "no answer".
53    pub lag1_autocorrelation: Option<T>,
54}
55
56/// Statistical test
57#[derive(Debug)]
58pub struct StatisticalTest<T: Float + Debug + Send + Sync + 'static> {
59    /// Test name
60    pub name: String,
61
62    /// Test type
63    pub test_type: StatisticalTestType,
64
65    /// Test statistic
66    pub test_statistic: T,
67
68    /// P-value
69    pub p_value: T,
70
71    /// Effect size
72    pub effect_size: T,
73
74    /// Confidence interval
75    pub confidence_interval: (T, T),
76}
77
78impl<T: Float + Debug + Default + std::iter::Sum + Send + Sync> StatisticalAnalyzer<T> {
79    pub(crate) fn new() -> Self {
80        Self {
81            statistical_tests: Vec::new(),
82            analysis_methods: vec![
83                AnalysisMethod::DescriptiveStatistics,
84                AnalysisMethod::CorrelationAnalysis,
85            ],
86            significance_thresholds: {
87                let mut thresholds = HashMap::new();
88                // Unreachable for `f32`/`f64` (see `crate::numeric`). If a
89                // caller's `Float` type did reject 0.05, a threshold of zero
90                // admits nothing — the conservative direction. A fallback of one
91                // would declare every test significant, which is precisely the
92                // silent-wrong-number failure `crate::numeric` exists to avoid.
93                thresholds.insert(
94                    "alpha".to_string(),
95                    scalar_as::<T>(0.05, "default significance level")
96                        .unwrap_or_else(|_| T::zero()),
97                );
98                thresholds
99            },
100            multiple_comparison: MultipleComparisonCorrection::BenjaminiHochberg,
101        }
102    }
103
104    /// Analyses this analyzer performs, in order.
105    pub fn analysis_methods(&self) -> &[AnalysisMethod] {
106        &self.analysis_methods
107    }
108
109    /// Choose which analyses [`StatisticalAnalyzer::analyze`] performs.
110    ///
111    /// Only the methods this analyzer actually implements are accepted. Enabling
112    /// one it cannot perform is reported as [`OptimError::NotImplemented`] rather
113    /// than accepted and silently skipped, which would let a caller believe a
114    /// regression or cluster analysis had run.
115    pub fn set_analysis_methods(&mut self, methods: Vec<AnalysisMethod>) -> Result<()> {
116        for method in &methods {
117            match method {
118                AnalysisMethod::DescriptiveStatistics | AnalysisMethod::CorrelationAnalysis => {}
119                other => {
120                    return Err(OptimError::NotImplemented(format!(
121                        "StatisticalAnalyzer implements DescriptiveStatistics and \
122                         CorrelationAnalysis; {:?} has no implementation here",
123                        other
124                    )))
125                }
126            }
127        }
128        self.analysis_methods = methods;
129        Ok(())
130    }
131
132    /// Run every enabled analysis over `values`.
133    pub fn analyze(&self, values: &[T]) -> Result<AnalysisReport<T>> {
134        let mut report = AnalysisReport {
135            descriptive: None,
136            lag1_autocorrelation: None,
137        };
138        for method in &self.analysis_methods {
139            match method {
140                AnalysisMethod::DescriptiveStatistics => {
141                    report.descriptive = Some(self.compute_descriptive_stats(values)?);
142                }
143                AnalysisMethod::CorrelationAnalysis => {
144                    report.lag1_autocorrelation = self.lag1_autocorrelation(values)?;
145                }
146                other => {
147                    return Err(OptimError::NotImplemented(format!(
148                        "StatisticalAnalyzer has no implementation for {:?}",
149                        other
150                    )))
151                }
152            }
153        }
154        Ok(report)
155    }
156
157    /// Lag-1 autocorrelation of `values`, or `None` when it is undefined.
158    fn lag1_autocorrelation(&self, values: &[T]) -> Result<Option<T>> {
159        if values.len() < 2 {
160            return Ok(None);
161        }
162        let count: T = count_as(values.len(), "sample size")?;
163        let mean = values.iter().cloned().sum::<T>() / count;
164
165        let mut denominator = T::zero();
166        for &value in values {
167            let centered = value - mean;
168            denominator = denominator + centered * centered;
169        }
170        if denominator <= T::zero() {
171            // A constant series has no autocorrelation to report.
172            return Ok(None);
173        }
174
175        let mut numerator = T::zero();
176        for window in values.windows(2) {
177            numerator = numerator + (window[0] - mean) * (window[1] - mean);
178        }
179        Ok(Some(numerator / denominator))
180    }
181
182    /// Recorded statistical tests, in the order they were added.
183    pub fn tests(&self) -> &[StatisticalTest<T>] {
184        &self.statistical_tests
185    }
186
187    /// Multiple-comparison correction applied by
188    /// [`StatisticalAnalyzer::adjusted_p_values`].
189    pub fn multiple_comparison(&self) -> &MultipleComparisonCorrection {
190        &self.multiple_comparison
191    }
192
193    /// Choose the multiple-comparison correction.
194    pub fn set_multiple_comparison(&mut self, correction: MultipleComparisonCorrection) {
195        self.multiple_comparison = correction;
196    }
197
198    /// p-values of the recorded tests after the configured multiple-comparison
199    /// correction, in the order the tests were added.
200    ///
201    /// Every correction is closed-form and clamped to `[0, 1]`:
202    ///
203    /// * `None` — the raw p-values.
204    /// * `Bonferroni` — `m * p`.
205    /// * `Sidak` — `1 - (1 - p)^m`.
206    /// * `HolmBonferroni` — step-down `(m - i) * p`, made monotone.
207    /// * `BenjaminiHochberg` — step-up `m * p / (i + 1)`, made monotone.
208    /// * `BenjaminiYekutieli` — Benjamini-Hochberg scaled by the harmonic number
209    ///   `sum_{k=1..m} 1/k`, which is the dependence-robust variant.
210    pub fn adjusted_p_values(&self) -> Result<Vec<T>> {
211        let m = self.statistical_tests.len();
212        if m == 0 {
213            return Ok(Vec::new());
214        }
215        let raw: Vec<T> = self
216            .statistical_tests
217            .iter()
218            .map(|test| test.p_value)
219            .collect();
220        let count: T = count_as(m, "number of statistical tests")?;
221
222        let adjusted = match self.multiple_comparison {
223            MultipleComparisonCorrection::None => raw,
224            MultipleComparisonCorrection::Bonferroni => {
225                raw.into_iter().map(|p| (p * count).min(T::one())).collect()
226            }
227            MultipleComparisonCorrection::Sidak => raw
228                .into_iter()
229                .map(|p| {
230                    let complement = (T::one() - p).max(T::zero());
231                    (T::one() - complement.powf(count)).min(T::one())
232                })
233                .collect(),
234            MultipleComparisonCorrection::HolmBonferroni => Self::step_down(&raw, |rank, p| {
235                let factor: T = count_as((m - rank).max(1), "Holm rank factor")?;
236                Ok((p * factor).min(T::one()))
237            })?,
238            MultipleComparisonCorrection::BenjaminiHochberg => {
239                Self::step_up(&raw, count, T::one())?
240            }
241            MultipleComparisonCorrection::BenjaminiYekutieli => {
242                let mut harmonic = T::zero();
243                for k in 1..=m {
244                    let denominator: T = count_as(k, "harmonic-series term")?;
245                    harmonic = harmonic + T::one() / denominator;
246                }
247                Self::step_up(&raw, count, harmonic)?
248            }
249        };
250        Ok(adjusted)
251    }
252
253    /// Indices of the tests whose adjusted p-value is at or below `alpha`.
254    pub fn significant_tests(&self, alpha: T) -> Result<Vec<usize>> {
255        let adjusted = self.adjusted_p_values()?;
256        Ok(adjusted
257            .into_iter()
258            .enumerate()
259            .filter(|(_, p)| *p <= alpha)
260            .map(|(index, _)| index)
261            .collect())
262    }
263
264    /// Shared machinery for the step-down (Holm) family: sort ascending, scale by
265    /// rank, then enforce monotonicity by carrying the running maximum forward.
266    fn step_down<F>(raw: &[T], scale: F) -> Result<Vec<T>>
267    where
268        F: Fn(usize, T) -> Result<T>,
269    {
270        let mut order: Vec<usize> = (0..raw.len()).collect();
271        order.sort_by(|&a, &b| raw[a].partial_cmp(&raw[b]).unwrap_or(Ordering::Equal));
272
273        let mut adjusted = vec![T::zero(); raw.len()];
274        let mut running = T::zero();
275        for (rank, &index) in order.iter().enumerate() {
276            let scaled = scale(rank, raw[index])?;
277            running = running.max(scaled);
278            adjusted[index] = running.min(T::one());
279        }
280        Ok(adjusted)
281    }
282
283    /// Shared machinery for the step-up (Benjamini) family: sort descending,
284    /// scale by `m * factor / rank`, then enforce monotonicity by carrying the
285    /// running minimum backwards.
286    fn step_up(raw: &[T], count: T, factor: T) -> Result<Vec<T>> {
287        let mut order: Vec<usize> = (0..raw.len()).collect();
288        order.sort_by(|&a, &b| raw[b].partial_cmp(&raw[a]).unwrap_or(Ordering::Equal));
289
290        let mut adjusted = vec![T::zero(); raw.len()];
291        let mut running = T::one();
292        for (position, &index) in order.iter().enumerate() {
293            let rank = raw.len() - position;
294            let rank_value: T = count_as(rank, "Benjamini rank")?;
295            let scaled = raw[index] * count * factor / rank_value;
296            running = running.min(scaled);
297            adjusted[index] = running.min(T::one()).max(T::zero());
298        }
299        Ok(adjusted)
300    }
301
302    pub(crate) fn compute_confidence_intervals(
303        &self,
304        results: &[TestResult<T>],
305    ) -> Result<HashMap<EvaluationMetric, (T, T)>> {
306        let mut intervals = HashMap::new();
307
308        if !results.is_empty() {
309            let scores: Vec<T> = results.iter().map(|r| r.score).collect();
310            let count: T = count_as(scores.len(), "confidence-interval sample size")?;
311            let mean = scores.iter().cloned().sum::<T>() / count;
312            let std_dev = if scores.len() > 1 {
313                let degrees_of_freedom: T =
314                    count_as(scores.len() - 1, "confidence-interval degrees of freedom")?;
315                let variance =
316                    scores.iter().map(|&s| (s - mean) * (s - mean)).sum::<T>() / degrees_of_freedom;
317                variance.sqrt()
318            } else {
319                T::zero()
320            };
321
322            // 95% confidence interval (normal approximation).
323            let z: T = scalar_as(1.96, "95% normal quantile")?;
324            let margin = std_dev * z / count.sqrt();
325            intervals.insert(
326                EvaluationMetric::FinalPerformance,
327                (mean - margin, mean + margin),
328            );
329        }
330
331        Ok(intervals)
332    }
333
334    /// Add a statistical test
335    pub fn add_test(&mut self, test: StatisticalTest<T>) {
336        self.statistical_tests.push(test);
337    }
338
339    /// Get significance threshold
340    pub fn get_threshold(&self, name: &str) -> Option<T> {
341        self.significance_thresholds.get(name).copied()
342    }
343
344    /// Set significance threshold
345    pub fn set_threshold(&mut self, name: String, value: T) {
346        self.significance_thresholds.insert(name, value);
347    }
348
349    /// Compute descriptive statistics.
350    ///
351    /// An empty sample yields an all-zero summary with `count == 0`, which the
352    /// caller can distinguish from a genuine all-zero sample. Fallible because the
353    /// sample size has to reach the element type to be a divisor; there is no
354    /// fallback divisor that would not silently report a different statistic.
355    pub fn compute_descriptive_stats(&self, values: &[T]) -> Result<DescriptiveStats<T>> {
356        if values.is_empty() {
357            return Ok(DescriptiveStats {
358                mean: T::zero(),
359                median: T::zero(),
360                std_dev: T::zero(),
361                min: T::zero(),
362                max: T::zero(),
363                count: 0,
364            });
365        }
366
367        let count = values.len();
368        let count_t: T = count_as(count, "descriptive sample size")?;
369        let mean = sum_of(values) / count_t;
370
371        let variance = values.iter().map(|&v| (v - mean) * (v - mean)).sum::<T>() / count_t;
372        let std_dev = variance.sqrt();
373
374        let mut sorted = values.to_vec();
375        sorted.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
376
377        let min = sorted.first().copied().unwrap_or(T::zero());
378        let max = sorted.last().copied().unwrap_or(T::zero());
379        let median = if count.is_multiple_of(2) {
380            (sorted[count / 2 - 1] + sorted[count / 2])
381                / scirs2_core::numeric::NumCast::from(2.0).unwrap_or_else(|| T::one())
382        } else {
383            sorted[count / 2]
384        };
385
386        Ok(DescriptiveStats {
387            mean,
388            median,
389            std_dev,
390            min,
391            max,
392            count,
393        })
394    }
395}
396
397/// Descriptive statistics result
398#[derive(Debug, Clone)]
399pub struct DescriptiveStats<T: Float + Debug + Send + Sync + 'static> {
400    /// Mean value
401    pub mean: T,
402    /// Median value
403    pub median: T,
404    /// Standard deviation
405    pub std_dev: T,
406    /// Minimum value
407    pub min: T,
408    /// Maximum value
409    pub max: T,
410    /// Count
411    pub count: usize,
412}
413
414/// Sum of a sample, factored out so the descriptive path reads the values once.
415fn sum_of<T: Float + std::iter::Sum>(values: &[T]) -> T {
416    values.iter().cloned().sum()
417}
418
419#[cfg(test)]
420mod tests {
421    use super::*;
422
423    fn test_with_p(name: &str, p_value: f64) -> StatisticalTest<f64> {
424        StatisticalTest {
425            name: name.to_string(),
426            test_type: StatisticalTestType::TTest,
427            test_statistic: 0.0,
428            p_value,
429            effect_size: 0.0,
430            confidence_interval: (0.0, 0.0),
431        }
432    }
433
434    fn analyzer_with(ps: &[f64]) -> StatisticalAnalyzer<f64> {
435        let mut analyzer = StatisticalAnalyzer::<f64>::new();
436        for (index, &p) in ps.iter().enumerate() {
437            analyzer.add_test(test_with_p(&format!("t{index}"), p));
438        }
439        analyzer
440    }
441
442    #[test]
443    fn no_correction_returns_the_raw_p_values() {
444        let mut analyzer = analyzer_with(&[0.01, 0.04, 0.20]);
445        analyzer.set_multiple_comparison(MultipleComparisonCorrection::None);
446        let adjusted = analyzer.adjusted_p_values().expect("adjustment");
447        assert_eq!(adjusted, vec![0.01, 0.04, 0.20]);
448    }
449
450    #[test]
451    fn bonferroni_and_sidak_multiply_by_the_family_size() {
452        let mut analyzer = analyzer_with(&[0.01, 0.04, 0.20]);
453
454        analyzer.set_multiple_comparison(MultipleComparisonCorrection::Bonferroni);
455        let bonferroni = analyzer.adjusted_p_values().expect("adjustment");
456        assert!((bonferroni[0] - 0.03).abs() < 1e-12);
457        assert!((bonferroni[1] - 0.12).abs() < 1e-12);
458        // Clamped at one rather than reporting a probability above 1.
459        assert!((bonferroni[2] - 0.6).abs() < 1e-12);
460
461        analyzer.set_multiple_comparison(MultipleComparisonCorrection::Sidak);
462        let sidak = analyzer.adjusted_p_values().expect("adjustment");
463        for (s, b) in sidak.iter().zip(bonferroni.iter()) {
464            // Sidak is uniformly less conservative than Bonferroni.
465            assert!(*s <= *b + 1e-12, "sidak {s} exceeded bonferroni {b}");
466            assert!((0.0..=1.0).contains(s));
467        }
468    }
469
470    #[test]
471    fn holm_is_monotone_and_less_conservative_than_bonferroni() {
472        let mut analyzer = analyzer_with(&[0.20, 0.01, 0.04]);
473        analyzer.set_multiple_comparison(MultipleComparisonCorrection::Bonferroni);
474        let bonferroni = analyzer.adjusted_p_values().expect("adjustment");
475
476        analyzer.set_multiple_comparison(MultipleComparisonCorrection::HolmBonferroni);
477        let holm = analyzer.adjusted_p_values().expect("adjustment");
478
479        // Smallest raw p (0.01, index 1) is scaled by 3, the next by 2, the last
480        // by 1 — then made monotone.
481        assert!((holm[1] - 0.03).abs() < 1e-12);
482        assert!((holm[2] - 0.08).abs() < 1e-12);
483        assert!((holm[0] - 0.20).abs() < 1e-12);
484        for (h, b) in holm.iter().zip(bonferroni.iter()) {
485            assert!(*h <= *b + 1e-12);
486        }
487    }
488
489    #[test]
490    fn benjamini_hochberg_is_monotone_and_yekutieli_is_stricter() {
491        let mut analyzer = analyzer_with(&[0.001, 0.008, 0.039, 0.041, 0.042]);
492
493        analyzer.set_multiple_comparison(MultipleComparisonCorrection::BenjaminiHochberg);
494        let bh = analyzer.adjusted_p_values().expect("adjustment");
495        // Ordering is preserved: a smaller raw p never gets a larger adjusted p.
496        for window in bh.windows(2) {
497            assert!(
498                window[0] <= window[1] + 1e-12,
499                "BH broke monotonicity: {bh:?}"
500            );
501        }
502        // 0.001 * 5 / 1 = 0.005.
503        assert!((bh[0] - 0.005).abs() < 1e-12);
504
505        analyzer.set_multiple_comparison(MultipleComparisonCorrection::BenjaminiYekutieli);
506        let by = analyzer.adjusted_p_values().expect("adjustment");
507        for (y, h) in by.iter().zip(bh.iter()) {
508            assert!(*y >= *h - 1e-12, "Yekutieli must not be looser than BH");
509            assert!((0.0..=1.0).contains(y));
510        }
511    }
512
513    #[test]
514    fn significance_uses_the_adjusted_values_not_the_raw_ones() {
515        let mut analyzer = analyzer_with(&[0.02, 0.03, 0.04]);
516        analyzer.set_multiple_comparison(MultipleComparisonCorrection::None);
517        assert_eq!(
518            analyzer.significant_tests(0.05).expect("significance"),
519            vec![0, 1, 2]
520        );
521
522        // Under Bonferroni every adjusted value is above 0.05, so nothing survives.
523        analyzer.set_multiple_comparison(MultipleComparisonCorrection::Bonferroni);
524        assert!(analyzer
525            .significant_tests(0.05)
526            .expect("significance")
527            .is_empty());
528
529        assert_eq!(analyzer.tests().len(), 3);
530        assert_eq!(
531            analyzer.multiple_comparison(),
532            &MultipleComparisonCorrection::Bonferroni
533        );
534    }
535
536    #[test]
537    fn an_empty_family_has_no_adjusted_values() {
538        let analyzer = StatisticalAnalyzer::<f64>::new();
539        assert!(analyzer.adjusted_p_values().expect("adjustment").is_empty());
540        assert!(analyzer
541            .significant_tests(0.05)
542            .expect("significance")
543            .is_empty());
544    }
545
546    #[test]
547    fn analyze_runs_exactly_the_enabled_methods() {
548        let mut analyzer = StatisticalAnalyzer::<f64>::new();
549        let series = [1.0, 2.0, 3.0, 4.0, 5.0];
550
551        let report = analyzer.analyze(&series).expect("analysis");
552        let descriptive = report.descriptive.expect("descriptive statistics enabled");
553        assert_eq!(descriptive.count, 5);
554        assert!((descriptive.mean - 3.0).abs() < 1e-12);
555        let autocorrelation = report
556            .lag1_autocorrelation
557            .expect("correlation analysis enabled");
558        // Standard (biased) lag-1 estimator on 1..=5: numerator 4, denominator 10.
559        assert!(
560            (autocorrelation - 0.4).abs() < 1e-12,
561            "got {autocorrelation}"
562        );
563
564        // An alternating series is negatively autocorrelated, so the sign is real
565        // signal and not an artefact of the ramp.
566        let alternating = [1.0, -1.0, 1.0, -1.0, 1.0, -1.0];
567        let flipped = analyzer
568            .analyze(&alternating)
569            .expect("analysis")
570            .lag1_autocorrelation
571            .expect("correlation analysis enabled");
572        assert!(flipped < 0.0, "got {flipped}");
573
574        analyzer
575            .set_analysis_methods(vec![AnalysisMethod::DescriptiveStatistics])
576            .expect("descriptive statistics are implemented");
577        let report = analyzer.analyze(&series).expect("analysis");
578        assert!(report.descriptive.is_some());
579        assert!(
580            report.lag1_autocorrelation.is_none(),
581            "a disabled analysis must not run"
582        );
583    }
584
585    #[test]
586    fn an_unimplemented_analysis_is_rejected_not_silently_skipped() {
587        let mut analyzer = StatisticalAnalyzer::<f64>::new();
588        let err = analyzer
589            .set_analysis_methods(vec![AnalysisMethod::PrincipalComponentAnalysis])
590            .expect_err("PCA has no implementation here");
591        assert!(matches!(err, OptimError::NotImplemented(_)));
592        // The rejected configuration must not have been stored.
593        assert_eq!(analyzer.analysis_methods().len(), 2);
594    }
595
596    #[test]
597    fn autocorrelation_is_undefined_for_short_or_constant_samples() {
598        let analyzer = StatisticalAnalyzer::<f64>::new();
599        assert!(analyzer
600            .analyze(&[])
601            .expect("analysis")
602            .lag1_autocorrelation
603            .is_none());
604        assert!(analyzer
605            .analyze(&[1.0])
606            .expect("analysis")
607            .lag1_autocorrelation
608            .is_none());
609        assert!(
610            analyzer
611                .analyze(&[2.0, 2.0, 2.0])
612                .expect("analysis")
613                .lag1_autocorrelation
614                .is_none(),
615            "a constant series has no autocorrelation to report"
616        );
617    }
618
619    #[test]
620    fn confidence_intervals_survive_degenerate_samples() {
621        let analyzer = StatisticalAnalyzer::<f64>::new();
622        // No results at all: no interval, and no panic from an empty divisor.
623        assert!(analyzer
624            .compute_confidence_intervals(&[])
625            .expect("intervals")
626            .is_empty());
627    }
628}