robust_rs/
estimator.rs

1//! The `RobustEstimator` trait (the sampling-theory result surface) and the
2//! shared regression fit type.
3
4use faer::linalg::solvers::DenseSolveCore;
5use ndarray::{Array1, Array2};
6use robust_rs_core::rho::RhoFunction;
7use robust_rs_core::types::Scale;
8
9const QUAD: usize = 128;
10
11/// A fitted robust regression, carrying enough state to report the theory.
12pub struct RegressionFit {
13    /// Estimated coefficients.
14    pub coefficients: Array1<f64>,
15    /// Estimated residual scale.
16    pub scale: Scale,
17    /// Residuals `y − Xβ̂`.
18    pub residuals: Array1<f64>,
19    /// Final IRLS weights.
20    pub weights: Array1<f64>,
21    /// The loss used, retained so influence/efficiency can be reported. This is
22    /// the estimator's **final-stage** loss (e.g. the M-step biweight for MM),
23    /// so the reported efficiency and covariance describe the delivered fit.
24    pub rho: Box<dyn RhoFunction>,
25    /// Asymptotic breakdown point ε*, carried as data set by the producing
26    /// estimator, *not* derived from `rho`, since ρ alone does not determine an
27    /// estimator's breakdown. Plain M-regression sets `0.0` (leverage); S/MM set
28    /// `min(δ/ρ_sup, 1 − δ/ρ_sup)`.
29    pub breakdown_point: f64,
30}
31
32/// Quantities every fitted robust estimator can report.
33pub trait RobustEstimator {
34    /// Estimated coefficients.
35    fn coefficients(&self) -> &Array1<f64>;
36    /// Estimated residual scale.
37    fn scale(&self) -> Scale;
38    /// The influence function `x ↦ ψ(x)/E[ψ']`.
39    fn influence_function(&self) -> Box<dyn Fn(f64) -> f64 + '_>;
40    /// Asymptotic variance `E[ψ²]/(E[ψ'])²`.
41    fn asymptotic_variance(&self) -> f64;
42    /// Efficiency relative to the Gaussian MLE.
43    fn gaussian_efficiency(&self) -> f64;
44    /// Approximate coefficient covariance `ŝ²·V·(XᵀX)⁻¹`.
45    fn coef_covariance(&self, x: &Array2<f64>) -> Array2<f64>;
46    /// Asymptotic breakdown point.
47    fn breakdown_point(&self) -> f64;
48}
49
50impl RobustEstimator for RegressionFit {
51    fn coefficients(&self) -> &Array1<f64> {
52        &self.coefficients
53    }
54    fn scale(&self) -> Scale {
55        self.scale
56    }
57    fn influence_function(&self) -> Box<dyn Fn(f64) -> f64 + '_> {
58        Box::new(robust_rs_core::theory::influence_function(&*self.rho, QUAD))
59    }
60    fn asymptotic_variance(&self) -> f64 {
61        robust_rs_core::theory::asymptotic_variance(&*self.rho, QUAD)
62    }
63    fn gaussian_efficiency(&self) -> f64 {
64        robust_rs_core::theory::gaussian_efficiency(&*self.rho, QUAD)
65    }
66    fn coef_covariance(&self, x: &Array2<f64>) -> Array2<f64> {
67        let p = x.ncols();
68        let s = self.scale.get();
69        let factor = s * s * self.asymptotic_variance();
70
71        // (XᵀX)⁻¹ via faer LU. XᵀX is SPD and full-rank after a successful fit.
72        let xtx = x.t().dot(x);
73        let a = faer::Mat::from_fn(p, p, |i, j| xtx[[i, j]]);
74        let inv = a.partial_piv_lu().inverse();
75
76        Array2::from_shape_fn((p, p), |(i, j)| inv[(i, j)] * factor)
77    }
78    fn breakdown_point(&self) -> f64 {
79        self.breakdown_point
80    }
81}
robust_rs/estimator.rs

robust_rs/
estimator.rs
