1use 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#[derive(Debug)]
26pub struct StatisticalAnalyzer<T: Float + Debug + Send + Sync + 'static> {
27 statistical_tests: Vec<StatisticalTest<T>>,
29
30 analysis_methods: Vec<AnalysisMethod>,
32
33 significance_thresholds: HashMap<String, T>,
35
36 multiple_comparison: MultipleComparisonCorrection,
39}
40
41#[derive(Debug, Clone)]
43pub struct AnalysisReport<T: Float + Debug + Send + Sync + 'static> {
44 pub descriptive: Option<DescriptiveStats<T>>,
47
48 pub lag1_autocorrelation: Option<T>,
54}
55
56#[derive(Debug)]
58pub struct StatisticalTest<T: Float + Debug + Send + Sync + 'static> {
59 pub name: String,
61
62 pub test_type: StatisticalTestType,
64
65 pub test_statistic: T,
67
68 pub p_value: T,
70
71 pub effect_size: T,
73
74 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 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 pub fn analysis_methods(&self) -> &[AnalysisMethod] {
106 &self.analysis_methods
107 }
108
109 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 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 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 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 pub fn tests(&self) -> &[StatisticalTest<T>] {
184 &self.statistical_tests
185 }
186
187 pub fn multiple_comparison(&self) -> &MultipleComparisonCorrection {
190 &self.multiple_comparison
191 }
192
193 pub fn set_multiple_comparison(&mut self, correction: MultipleComparisonCorrection) {
195 self.multiple_comparison = correction;
196 }
197
198 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 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 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 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 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 pub fn add_test(&mut self, test: StatisticalTest<T>) {
336 self.statistical_tests.push(test);
337 }
338
339 pub fn get_threshold(&self, name: &str) -> Option<T> {
341 self.significance_thresholds.get(name).copied()
342 }
343
344 pub fn set_threshold(&mut self, name: String, value: T) {
346 self.significance_thresholds.insert(name, value);
347 }
348
349 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#[derive(Debug, Clone)]
399pub struct DescriptiveStats<T: Float + Debug + Send + Sync + 'static> {
400 pub mean: T,
402 pub median: T,
404 pub std_dev: T,
406 pub min: T,
408 pub max: T,
410 pub count: usize,
412}
413
414fn 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 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 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 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 for window in bh.windows(2) {
497 assert!(
498 window[0] <= window[1] + 1e-12,
499 "BH broke monotonicity: {bh:?}"
500 );
501 }
502 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 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 assert!(
560 (autocorrelation - 0.4).abs() < 1e-12,
561 "got {autocorrelation}"
562 );
563
564 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 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 assert!(analyzer
624 .compute_confidence_intervals(&[])
625 .expect("intervals")
626 .is_empty());
627 }
628}