1use scirs2_core::ndarray::{s, Array1, Array2, Axis};
8use scirs2_core::random::{prelude::*, SeedableRng};
9use sklears_core::{
10 error::{Result, SklearsError},
11 types::Float,
12};
13
14#[derive(Debug, Clone)]
16pub struct UncertaintyConfig {
17 pub n_bootstrap_samples: usize,
19 pub n_mc_samples: usize,
21 pub confidence_level: Float,
23 pub method: UncertaintyMethod,
25 pub decompose_uncertainty: bool,
27 pub random_seed: Option<u64>,
29}
30
31impl Default for UncertaintyConfig {
32 fn default() -> Self {
33 Self {
34 n_bootstrap_samples: 1000,
35 n_mc_samples: 100,
36 confidence_level: 0.95,
37 method: UncertaintyMethod::Bootstrap,
38 decompose_uncertainty: true,
39 random_seed: None,
40 }
41 }
42}
43
44#[derive(Debug, Clone)]
46pub enum UncertaintyMethod {
47 Bootstrap,
49 Bayesian,
51 Conformal,
53 Ensemble { n_models: usize },
55 Dropout {
57 dropout_rate: Float,
58 n_forward_passes: usize,
59 },
60}
61
62#[derive(Debug, Clone)]
64pub struct UncertaintyResult {
65 pub predictions: Array1<Float>,
67 pub total_uncertainty: Array1<Float>,
69 pub epistemic_uncertainty: Option<Array1<Float>>,
71 pub aleatoric_uncertainty: Option<Array1<Float>>,
73 pub lower_bound: Array1<Float>,
75 pub upper_bound: Array1<Float>,
77 pub confidence_level: Float,
79 pub method: String,
81}
82
83impl UncertaintyResult {
84 pub fn interval_width(&self) -> Array1<Float> {
86 &self.upper_bound - &self.lower_bound
87 }
88
89 pub fn coverage(&self, y_true: &Array1<Float>) -> Result<Float> {
91 if y_true.len() != self.predictions.len() {
92 return Err(SklearsError::DimensionMismatch {
93 expected: self.predictions.len(),
94 actual: y_true.len(),
95 });
96 }
97
98 let in_interval = y_true
99 .iter()
100 .zip(self.lower_bound.iter())
101 .zip(self.upper_bound.iter())
102 .map(|((&y, &lower), &upper)| if y >= lower && y <= upper { 1.0 } else { 0.0 })
103 .sum::<Float>();
104
105 Ok(in_interval / y_true.len() as Float)
106 }
107
108 pub fn mean_interval_width(&self) -> Float {
110 self.interval_width().mean().unwrap_or(0.0)
111 }
112}
113
114#[derive(Debug)]
116pub struct UncertaintyQuantifier {
117 config: UncertaintyConfig,
118}
119
120impl UncertaintyQuantifier {
121 pub fn new(config: UncertaintyConfig) -> Self {
123 Self { config }
124 }
125
126 pub fn bootstrap_uncertainty<F>(
128 &self,
129 X_train: &Array2<Float>,
130 y_train: &Array1<Float>,
131 X_test: &Array2<Float>,
132 fit_predict_fn: F,
133 ) -> Result<UncertaintyResult>
134 where
135 F: Fn(&Array2<Float>, &Array1<Float>, &Array2<Float>) -> Result<Array1<Float>>,
136 {
137 let mut rng = self.create_rng();
138 let n_train = X_train.nrows();
139 let n_test = X_test.nrows();
140
141 let mut bootstrap_predictions = Array2::zeros((self.config.n_bootstrap_samples, n_test));
143
144 for i in 0..self.config.n_bootstrap_samples {
145 let bootstrap_indices = self.bootstrap_sample(n_train, &mut rng);
147 let (X_boot, y_boot) =
148 self.extract_bootstrap_data(X_train, y_train, &bootstrap_indices);
149
150 let predictions = fit_predict_fn(&X_boot, &y_boot, X_test)?;
152 bootstrap_predictions
153 .slice_mut(s![i, ..])
154 .assign(&predictions);
155 }
156
157 self.compute_bootstrap_statistics(&bootstrap_predictions)
159 }
160
161 pub fn bayesian_uncertainty(
163 &self,
164 posterior_samples: &Array2<Float>,
165 X_test: &Array2<Float>,
166 noise_precision: Float,
167 ) -> Result<UncertaintyResult> {
168 let n_samples = posterior_samples.nrows();
169 let n_test = X_test.nrows();
170
171 let mut predictions = Array2::zeros((n_samples, n_test));
173
174 for i in 0..n_samples {
175 let weights = posterior_samples.slice(s![i, ..]);
176 let pred = X_test.dot(&weights);
177 predictions.slice_mut(s![i, ..]).assign(&pred);
178 }
179
180 let pred_mean = predictions.mean_axis(Axis(0)).ok_or_else(|| {
182 SklearsError::NumericalError(
183 "mean computation should succeed for non-empty array".into(),
184 )
185 })?;
186 let pred_var = predictions.var_axis(Axis(0), 0.0);
187
188 let epistemic_var = pred_var.clone(); let aleatoric_var = Array1::from_elem(n_test, 1.0 / noise_precision); let total_var = &epistemic_var + &aleatoric_var;
192
193 let total_std = total_var.mapv(|v| v.sqrt());
194 let epistemic_std = epistemic_var.mapv(|v| v.sqrt());
195 let aleatoric_std = aleatoric_var.mapv(|v| v.sqrt());
196
197 let alpha = 1.0 - self.config.confidence_level;
199 let z_score = self.compute_normal_quantile(1.0 - alpha / 2.0)?;
200
201 let lower_bound = &pred_mean - z_score * &total_std;
202 let upper_bound = &pred_mean + z_score * &total_std;
203
204 Ok(UncertaintyResult {
205 predictions: pred_mean,
206 total_uncertainty: total_std,
207 epistemic_uncertainty: Some(epistemic_std),
208 aleatoric_uncertainty: Some(aleatoric_std),
209 lower_bound,
210 upper_bound,
211 confidence_level: self.config.confidence_level,
212 method: "Bayesian".to_string(),
213 })
214 }
215
216 pub fn conformal_uncertainty<F>(
218 &self,
219 X_cal: &Array2<Float>,
220 y_cal: &Array1<Float>,
221 X_test: &Array2<Float>,
222 predict_fn: F,
223 ) -> Result<UncertaintyResult>
224 where
225 F: Fn(&Array2<Float>) -> Result<Array1<Float>>,
226 {
227 let cal_predictions = predict_fn(X_cal)?;
229 let cal_residuals: Array1<Float> = (y_cal - &cal_predictions).mapv(|r| r.abs());
230
231 let alpha = 1.0 - self.config.confidence_level;
233 let quantile_level = 1.0 - alpha;
234 let quantile = self.compute_empirical_quantile(&cal_residuals, quantile_level)?;
235
236 let test_predictions = predict_fn(X_test)?;
238
239 let lower_bound = &test_predictions - quantile;
241 let upper_bound = &test_predictions + quantile;
242
243 let total_uncertainty = Array1::from_elem(test_predictions.len(), quantile);
245
246 Ok(UncertaintyResult {
247 predictions: test_predictions,
248 total_uncertainty,
249 epistemic_uncertainty: None,
250 aleatoric_uncertainty: None,
251 lower_bound,
252 upper_bound,
253 confidence_level: self.config.confidence_level,
254 method: "Conformal".to_string(),
255 })
256 }
257
258 pub fn ensemble_uncertainty(
260 &self,
261 ensemble_predictions: &Array2<Float>,
262 ) -> Result<UncertaintyResult> {
263 let pred_mean = ensemble_predictions.mean_axis(Axis(0)).ok_or_else(|| {
264 SklearsError::NumericalError(
265 "mean computation should succeed for non-empty array".into(),
266 )
267 })?;
268 let pred_var = ensemble_predictions.var_axis(Axis(0), 0.0);
269 let pred_std = pred_var.mapv(|v| v.sqrt());
270
271 let df = ensemble_predictions.nrows() - 1;
273 let alpha = 1.0 - self.config.confidence_level;
274 let t_critical = self.compute_t_quantile(1.0 - alpha / 2.0, df as Float)?;
275
276 let margin = t_critical * &pred_std / (ensemble_predictions.nrows() as Float).sqrt();
277 let lower_bound = &pred_mean - &margin;
278 let upper_bound = &pred_mean + &margin;
279
280 Ok(UncertaintyResult {
281 predictions: pred_mean,
282 total_uncertainty: pred_std.clone(),
283 epistemic_uncertainty: Some(pred_std), aleatoric_uncertainty: None,
285 lower_bound,
286 upper_bound,
287 confidence_level: self.config.confidence_level,
288 method: "Ensemble".to_string(),
289 })
290 }
291
292 fn create_rng(&self) -> StdRng {
294 if let Some(seed) = self.config.random_seed {
295 StdRng::seed_from_u64(seed)
296 } else {
297 StdRng::from_rng(&mut scirs2_core::random::thread_rng())
298 }
299 }
300
301 fn bootstrap_sample(&self, n_samples: usize, rng: &mut impl Rng) -> Vec<usize> {
303 (0..n_samples)
304 .map(|_| rng.random_range(0..n_samples))
305 .collect()
306 }
307
308 fn extract_bootstrap_data(
310 &self,
311 X: &Array2<Float>,
312 y: &Array1<Float>,
313 indices: &[usize],
314 ) -> (Array2<Float>, Array1<Float>) {
315 let n_samples = indices.len();
316 let n_features = X.ncols();
317
318 let mut X_boot = Array2::zeros((n_samples, n_features));
319 let mut y_boot = Array1::zeros(n_samples);
320
321 for (i, &idx) in indices.iter().enumerate() {
322 X_boot.slice_mut(s![i, ..]).assign(&X.slice(s![idx, ..]));
323 y_boot[i] = y[idx];
324 }
325
326 (X_boot, y_boot)
327 }
328
329 fn compute_bootstrap_statistics(
331 &self,
332 predictions: &Array2<Float>,
333 ) -> Result<UncertaintyResult> {
334 let pred_mean = predictions.mean_axis(Axis(0)).ok_or_else(|| {
335 SklearsError::NumericalError(
336 "mean computation should succeed for non-empty array".into(),
337 )
338 })?;
339 let pred_std = predictions.std_axis(Axis(0), 0.0);
340
341 let alpha = 1.0 - self.config.confidence_level;
343 let lower_quantile = alpha / 2.0;
344 let upper_quantile = 1.0 - alpha / 2.0;
345
346 let n_test = predictions.ncols();
347 let mut lower_bound = Array1::zeros(n_test);
348 let mut upper_bound = Array1::zeros(n_test);
349
350 for i in 0..n_test {
351 let col = predictions.column(i);
352 let mut sorted_col: Vec<Float> = col.to_vec();
353 sorted_col.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
354
355 lower_bound[i] =
356 self.compute_empirical_quantile_from_sorted(&sorted_col, lower_quantile)?;
357 upper_bound[i] =
358 self.compute_empirical_quantile_from_sorted(&sorted_col, upper_quantile)?;
359 }
360
361 Ok(UncertaintyResult {
362 predictions: pred_mean,
363 total_uncertainty: pred_std.clone(),
364 epistemic_uncertainty: Some(pred_std), aleatoric_uncertainty: None,
366 lower_bound,
367 upper_bound,
368 confidence_level: self.config.confidence_level,
369 method: "Bootstrap".to_string(),
370 })
371 }
372
373 fn compute_empirical_quantile(&self, data: &Array1<Float>, quantile: Float) -> Result<Float> {
375 if !(0.0..=1.0).contains(&quantile) {
376 return Err(SklearsError::InvalidParameter {
377 name: "quantile".to_string(),
378 reason: format!("Quantile must be between 0 and 1, got {}", quantile),
379 });
380 }
381
382 let mut sorted_data: Vec<Float> = data.to_vec();
383 sorted_data.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
384
385 self.compute_empirical_quantile_from_sorted(&sorted_data, quantile)
386 }
387
388 fn compute_empirical_quantile_from_sorted(
390 &self,
391 sorted_data: &[Float],
392 quantile: Float,
393 ) -> Result<Float> {
394 if sorted_data.is_empty() {
395 return Err(SklearsError::InvalidInput(
396 "Empty data for quantile computation".to_string(),
397 ));
398 }
399
400 let n = sorted_data.len();
401 let index = quantile * (n - 1) as Float;
402 let lower_index = index.floor() as usize;
403 let upper_index = index.ceil() as usize;
404
405 if lower_index == upper_index {
406 Ok(sorted_data[lower_index])
407 } else {
408 let weight = index - lower_index as Float;
409 Ok(sorted_data[lower_index] * (1.0 - weight) + sorted_data[upper_index] * weight)
410 }
411 }
412
413 #[allow(clippy::only_used_in_recursion)]
415 fn compute_normal_quantile(&self, p: Float) -> Result<Float> {
416 if p <= 0.0 || p >= 1.0 {
417 return Err(SklearsError::InvalidParameter {
418 name: "p".to_string(),
419 reason: format!("Probability must be between 0 and 1, got {}", p),
420 });
421 }
422
423 let a = [
426 -3.969683028665376e+01,
427 2.209460984245205e+02,
428 -2.759285104469687e+02,
429 1.383_577_518_672_69e2,
430 -3.066479806614716e+01,
431 2.506628277459239e+00,
432 ];
433 let b = [
434 -5.447609879822406e+01,
435 1.615858368580409e+02,
436 -1.556989798598866e+02,
437 6.680131188771972e+01,
438 -1.328068155288572e+01,
439 1.0,
440 ];
441
442 let result = if p == 0.5 {
443 0.0
444 } else if p > 0.5 {
445 let q = p - 0.5;
447 let r = q * q;
448 let num = a[5] + r * (a[4] + r * (a[3] + r * (a[2] + r * (a[1] + r * a[0]))));
449 let den = b[5] + r * (b[4] + r * (b[3] + r * (b[2] + r * (b[1] + r * b[0]))));
450 q * num / den
451 } else {
452 -self.compute_normal_quantile(1.0 - p)?
454 };
455
456 Ok(result)
457 }
458
459 fn compute_t_quantile(&self, p: Float, df: Float) -> Result<Float> {
461 if df <= 0.0 {
462 return Err(SklearsError::InvalidParameter {
463 name: "df".to_string(),
464 reason: format!("Degrees of freedom must be positive, got {}", df),
465 });
466 }
467
468 if df >= 30.0 {
470 return self.compute_normal_quantile(p);
471 }
472
473 let z = self.compute_normal_quantile(p)?;
475 let correction = z * z * z / (4.0 * df) + z * z * z * z * z / (96.0 * df * df);
476
477 Ok(z + correction)
478 }
479}
480
481impl Default for UncertaintyQuantifier {
482 fn default() -> Self {
483 Self::new(UncertaintyConfig::default())
484 }
485}
486
487pub trait UncertaintyCapable {
489 fn predict_with_uncertainty(
491 &self,
492 X: &Array2<Float>,
493 config: &UncertaintyConfig,
494 ) -> Result<UncertaintyResult>;
495
496 fn epistemic_uncertainty(&self, X: &Array2<Float>) -> Result<Array1<Float>>;
498
499 fn aleatoric_uncertainty(&self, X: &Array2<Float>) -> Result<Array1<Float>>;
501}
502
503#[derive(Debug, Clone)]
505pub struct CalibrationMetrics {
506 pub expected_calibration_error: Float,
508 pub maximum_calibration_error: Float,
510 pub brier_score: Option<Float>,
512 pub coverage: Float,
514 pub mean_interval_width: Float,
516}
517
518impl CalibrationMetrics {
519 pub fn compute(
521 uncertainty_result: &UncertaintyResult,
522 y_true: &Array1<Float>,
523 _n_bins: usize,
524 ) -> Result<Self> {
525 let coverage = uncertainty_result.coverage(y_true)?;
526 let mean_interval_width = uncertainty_result.mean_interval_width();
527
528 let expected_calibration_error = (coverage - uncertainty_result.confidence_level).abs();
530 let maximum_calibration_error = expected_calibration_error; Ok(Self {
533 expected_calibration_error,
534 maximum_calibration_error,
535 brier_score: None, coverage,
537 mean_interval_width,
538 })
539 }
540}
541
542#[allow(non_snake_case)]
543#[cfg(test)]
544mod tests {
545 use super::*;
546 use scirs2_core::ndarray::Array;
547
548 #[test]
549 fn test_uncertainty_config() {
550 let config = UncertaintyConfig::default();
551 assert_eq!(config.n_bootstrap_samples, 1000);
552 assert_eq!(config.confidence_level, 0.95);
553 assert!(config.decompose_uncertainty);
554 }
555
556 #[test]
557 fn test_uncertainty_result() {
558 let predictions = Array::from_vec(vec![1.0, 2.0, 3.0]);
559 let lower_bound = Array::from_vec(vec![0.5, 1.5, 2.5]);
560 let upper_bound = Array::from_vec(vec![1.5, 2.5, 3.5]);
561
562 let result = UncertaintyResult {
563 predictions: predictions.clone(),
564 total_uncertainty: Array::from_vec(vec![0.25, 0.25, 0.25]),
565 epistemic_uncertainty: None,
566 aleatoric_uncertainty: None,
567 lower_bound: lower_bound.clone(),
568 upper_bound: upper_bound.clone(),
569 confidence_level: 0.95,
570 method: "Test".to_string(),
571 };
572
573 let width = result.interval_width();
574 assert_eq!(width[0], 1.0);
575 assert_eq!(width[1], 1.0);
576 assert_eq!(width[2], 1.0);
577
578 assert_eq!(result.mean_interval_width(), 1.0);
579
580 let y_true = Array::from_vec(vec![1.2, 2.1, 2.8]);
582 let coverage = result.coverage(&y_true).expect("operation should succeed");
583 assert_eq!(coverage, 1.0); }
585
586 #[test]
587 fn test_empirical_quantile() {
588 let quantifier = UncertaintyQuantifier::default();
589 let data = Array::from_vec(vec![1.0, 2.0, 3.0, 4.0, 5.0]);
590
591 let median = quantifier
592 .compute_empirical_quantile(&data, 0.5)
593 .expect("operation should succeed");
594 assert_eq!(median, 3.0);
595
596 let q25 = quantifier
597 .compute_empirical_quantile(&data, 0.25)
598 .expect("operation should succeed");
599 assert_eq!(q25, 2.0);
600
601 let q75 = quantifier
602 .compute_empirical_quantile(&data, 0.75)
603 .expect("operation should succeed");
604 assert_eq!(q75, 4.0);
605 }
606
607 #[test]
608 fn test_normal_quantile() {
609 let quantifier = UncertaintyQuantifier::default();
610
611 let median = quantifier
613 .compute_normal_quantile(0.5)
614 .expect("operation should succeed");
615 assert!(median.abs() < 1e-10);
616
617 let q975 = quantifier
619 .compute_normal_quantile(0.975)
620 .expect("operation should succeed");
621 assert!((q975 - 1.96).abs() < 0.2, "Expected ~1.96, got {}", q975);
623 }
624
625 #[test]
626 fn test_bootstrap_sample() {
627 let quantifier = UncertaintyQuantifier::new(UncertaintyConfig {
628 random_seed: Some(42),
629 ..Default::default()
630 });
631 let mut rng = quantifier.create_rng();
632
633 let indices = quantifier.bootstrap_sample(5, &mut rng);
634 assert_eq!(indices.len(), 5);
635 assert!(indices.iter().all(|&i| i < 5));
636 }
637
638 #[test]
639 #[allow(non_snake_case)]
640 fn test_bayesian_uncertainty() {
641 let quantifier = UncertaintyQuantifier::default();
642
643 let posterior_samples = Array::from_shape_vec((3, 2), vec![1.0, 0.5, 1.1, 0.4, 0.9, 0.6])
645 .expect("valid array shape");
646
647 let X_test =
649 Array::from_shape_vec((2, 2), vec![1.0, 1.0, 2.0, 2.0]).expect("valid array shape");
650
651 let noise_precision = 1.0;
652
653 let result = quantifier
654 .bayesian_uncertainty(&posterior_samples, &X_test, noise_precision)
655 .expect("operation should succeed");
656
657 assert_eq!(result.predictions.len(), 2);
658 assert_eq!(result.total_uncertainty.len(), 2);
659 assert!(result.epistemic_uncertainty.is_some());
660 assert!(result.aleatoric_uncertainty.is_some());
661 assert_eq!(result.method, "Bayesian");
662 }
663
664 #[test]
665 fn test_ensemble_uncertainty() {
666 let quantifier = UncertaintyQuantifier::default();
667
668 let ensemble_predictions = Array::from_shape_vec(
670 (5, 3),
671 vec![
672 1.0, 2.0, 3.0, 1.1, 1.9, 3.1, 0.9, 2.1, 2.9, 1.0, 2.0, 3.0, 1.0, 2.0, 3.0,
673 ],
674 )
675 .expect("operation should succeed");
676
677 let result = quantifier
678 .ensemble_uncertainty(&ensemble_predictions)
679 .expect("operation should succeed");
680
681 assert_eq!(result.predictions.len(), 3);
682 assert_eq!(result.total_uncertainty.len(), 3);
683 assert!(result.epistemic_uncertainty.is_some());
684 assert_eq!(result.method, "Ensemble");
685 }
686
687 #[test]
688 fn test_calibration_metrics() {
689 let uncertainty_result = UncertaintyResult {
690 predictions: Array::from_vec(vec![1.0, 2.0, 3.0]),
691 total_uncertainty: Array::from_vec(vec![0.1, 0.1, 0.1]),
692 epistemic_uncertainty: None,
693 aleatoric_uncertainty: None,
694 lower_bound: Array::from_vec(vec![0.8, 1.8, 2.8]),
695 upper_bound: Array::from_vec(vec![1.2, 2.2, 3.2]),
696 confidence_level: 0.95,
697 method: "Test".to_string(),
698 };
699
700 let y_true = Array::from_vec(vec![1.1, 2.1, 3.1]);
701
702 let metrics = CalibrationMetrics::compute(&uncertainty_result, &y_true, 10)
703 .expect("operation should succeed");
704
705 assert_eq!(metrics.coverage, 1.0);
706 assert!((metrics.mean_interval_width - 0.4).abs() < 1e-10);
707 assert!((metrics.expected_calibration_error - 0.05).abs() < 1e-10);
708 }
709}