Skip to main content

sklears_linear/
uncertainty_quantification.rs

1//! Uncertainty Quantification for Linear Models
2//!
3//! This module provides comprehensive uncertainty quantification methods for linear models,
4//! including epistemic and aleatoric uncertainty, conformal prediction, and Bayesian
5//! uncertainty propagation.
6
7use 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/// Configuration for uncertainty quantification
15#[derive(Debug, Clone)]
16pub struct UncertaintyConfig {
17    /// Number of bootstrap samples for uncertainty estimation
18    pub n_bootstrap_samples: usize,
19    /// Number of Monte Carlo samples for Bayesian uncertainty
20    pub n_mc_samples: usize,
21    /// Confidence level for prediction intervals (e.g., 0.95 for 95% intervals)
22    pub confidence_level: Float,
23    /// Method for uncertainty quantification
24    pub method: UncertaintyMethod,
25    /// Whether to separate epistemic and aleatoric uncertainty
26    pub decompose_uncertainty: bool,
27    /// Random seed for reproducibility
28    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/// Methods for uncertainty quantification
45#[derive(Debug, Clone)]
46pub enum UncertaintyMethod {
47    /// Bootstrap-based uncertainty quantification
48    Bootstrap,
49    /// Bayesian uncertainty propagation
50    Bayesian,
51    /// Conformal prediction
52    Conformal,
53    /// Ensemble-based uncertainty
54    Ensemble { n_models: usize },
55    /// Dropout-based uncertainty (for neural approaches)
56    Dropout {
57        dropout_rate: Float,
58        n_forward_passes: usize,
59    },
60}
61
62/// Results of uncertainty quantification
63#[derive(Debug, Clone)]
64pub struct UncertaintyResult {
65    /// Predicted mean values
66    pub predictions: Array1<Float>,
67    /// Total predictive uncertainty (standard deviation)
68    pub total_uncertainty: Array1<Float>,
69    /// Epistemic uncertainty (model uncertainty)
70    pub epistemic_uncertainty: Option<Array1<Float>>,
71    /// Aleatoric uncertainty (data uncertainty)
72    pub aleatoric_uncertainty: Option<Array1<Float>>,
73    /// Lower bound of prediction interval
74    pub lower_bound: Array1<Float>,
75    /// Upper bound of prediction interval
76    pub upper_bound: Array1<Float>,
77    /// Confidence level used
78    pub confidence_level: Float,
79    /// Method used for uncertainty quantification
80    pub method: String,
81}
82
83impl UncertaintyResult {
84    /// Compute prediction interval width
85    pub fn interval_width(&self) -> Array1<Float> {
86        &self.upper_bound - &self.lower_bound
87    }
88
89    /// Check if true values fall within prediction intervals
90    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    /// Compute mean interval width
109    pub fn mean_interval_width(&self) -> Float {
110        self.interval_width().mean().unwrap_or(0.0)
111    }
112}
113
114/// Uncertainty quantification engine
115#[derive(Debug)]
116pub struct UncertaintyQuantifier {
117    config: UncertaintyConfig,
118}
119
120impl UncertaintyQuantifier {
121    /// Create a new uncertainty quantifier
122    pub fn new(config: UncertaintyConfig) -> Self {
123        Self { config }
124    }
125
126    /// Quantify uncertainty using bootstrap method
127    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        // Collect bootstrap predictions
142        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            // Create bootstrap sample
146            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            // Fit model and predict
151            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        // Compute statistics
158        self.compute_bootstrap_statistics(&bootstrap_predictions)
159    }
160
161    /// Quantify uncertainty using Bayesian approach
162    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        // Compute predictions for each posterior sample
172        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        // Compute predictive statistics
181        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        // Decompose uncertainty
189        let epistemic_var = pred_var.clone(); // Model uncertainty
190        let aleatoric_var = Array1::from_elem(n_test, 1.0 / noise_precision); // Noise uncertainty
191        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        // Compute prediction intervals
198        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    /// Quantify uncertainty using conformal prediction
217    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        // Compute calibration residuals
228        let cal_predictions = predict_fn(X_cal)?;
229        let cal_residuals: Array1<Float> = (y_cal - &cal_predictions).mapv(|r| r.abs());
230
231        // Compute conformal quantile
232        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        // Make predictions on test set
237        let test_predictions = predict_fn(X_test)?;
238
239        // Construct prediction intervals
240        let lower_bound = &test_predictions - quantile;
241        let upper_bound = &test_predictions + quantile;
242
243        // Conformal prediction provides marginal coverage, not pointwise uncertainty
244        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    /// Quantify uncertainty using ensemble approach
259    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        // Compute prediction intervals using t-distribution for small ensembles
272        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), // Ensemble captures epistemic uncertainty
284            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    /// Create random number generator
293    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    /// Generate bootstrap sample indices
302    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    /// Extract bootstrap data based on indices
309    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    /// Compute statistics from bootstrap predictions
330    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        // Compute empirical quantiles for prediction intervals
342        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), // Bootstrap captures epistemic uncertainty
365            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    /// Compute empirical quantile from data
374    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    /// Compute empirical quantile from sorted data
389    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    /// Compute normal distribution quantile (approximate)
414    #[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        // Use a simpler and more reliable approximation
424        // Acklam's approximation for the normal quantile function
425        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            // Upper half
446            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            // Lower half - use symmetry
453            -self.compute_normal_quantile(1.0 - p)?
454        };
455
456        Ok(result)
457    }
458
459    /// Compute t-distribution quantile (approximate)
460    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        // For large df, t-distribution approaches normal
469        if df >= 30.0 {
470            return self.compute_normal_quantile(p);
471        }
472
473        // Simple approximation for t-quantile
474        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
487/// Trait for models that support uncertainty quantification
488pub trait UncertaintyCapable {
489    /// Predict with uncertainty quantification
490    fn predict_with_uncertainty(
491        &self,
492        X: &Array2<Float>,
493        config: &UncertaintyConfig,
494    ) -> Result<UncertaintyResult>;
495
496    /// Get epistemic uncertainty (model uncertainty)
497    fn epistemic_uncertainty(&self, X: &Array2<Float>) -> Result<Array1<Float>>;
498
499    /// Get aleatoric uncertainty (data uncertainty)
500    fn aleatoric_uncertainty(&self, X: &Array2<Float>) -> Result<Array1<Float>>;
501}
502
503/// Calibration metrics for uncertainty quantification
504#[derive(Debug, Clone)]
505pub struct CalibrationMetrics {
506    /// Expected calibration error
507    pub expected_calibration_error: Float,
508    /// Maximum calibration error
509    pub maximum_calibration_error: Float,
510    /// Brier score (for probabilistic predictions)
511    pub brier_score: Option<Float>,
512    /// Coverage probability
513    pub coverage: Float,
514    /// Mean interval width
515    pub mean_interval_width: Float,
516}
517
518impl CalibrationMetrics {
519    /// Compute calibration metrics for uncertainty predictions
520    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        // Compute calibration error (simplified implementation)
529        let expected_calibration_error = (coverage - uncertainty_result.confidence_level).abs();
530        let maximum_calibration_error = expected_calibration_error; // Simplified
531
532        Ok(Self {
533            expected_calibration_error,
534            maximum_calibration_error,
535            brier_score: None, // Would need probabilistic predictions
536            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        // Test coverage
581        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); // All points within intervals
584    }
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        // Test median (should be approximately 0)
612        let median = quantifier
613            .compute_normal_quantile(0.5)
614            .expect("operation should succeed");
615        assert!(median.abs() < 1e-10);
616
617        // Test 97.5% quantile (should be approximately 1.96)
618        let q975 = quantifier
619            .compute_normal_quantile(0.975)
620            .expect("operation should succeed");
621        // More relaxed tolerance for numerical approximation algorithms
622        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        // Mock posterior samples (3 samples, 2 features)
644        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        // Test data (2 samples, 2 features)
648        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        // Mock ensemble predictions (5 models, 3 test points)
669        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}