robust_rs_core/scale/
s_scale.rs

1//! The S-estimator of scale.
2
3use crate::error::RobustError;
4use crate::rho::RhoFunction;
5use crate::scale::ScaleEstimator;
6use crate::solver::Control;
7use crate::theory::gauss_hermite;
8use crate::types::Scale;
9
10/// The M-estimate of scale ("S-scale"): the `s > 0` solving
11/// `(1/n) Σ ρ(rᵢ/s) = δ`.
12///
13/// With a bounded, redescending `ρ` and `δ = ρ.rho_sup()/2` this is the
14/// 50%-breakdown scale that S- and MM-regression minimize; with `δ = E_Φ[ρ]` it
15/// is Fisher-consistent for σ at the Gaussian (for the biweight the S-step
16/// tuning `c ≈ 1.547` makes the two coincide). Generic over the loss, but the
17/// loss **must** have bounded `ρ` (`rho_sup().is_some()`): an unbounded `ρ`
18/// gives zero breakdown (a single `rᵢ → ∞` drives `s` without limit) so
19/// `new` rejects it with `RobustError::UnboundedLoss`.
20///
21/// Solved by the standard fixed-point iteration
22/// `s ← s·√(mean_i ρ(rᵢ/s) / δ)` from a MAD start; `g(s) = (1/n)Σρ(rᵢ/s)` is
23/// monotone decreasing, so the root is unique.
24#[derive(Debug, Clone, Copy)]
25pub struct SScale<R: RhoFunction> {
26    rho: R,
27    delta: f64,
28    control: Control,
29}
30
31impl<R: RhoFunction> SScale<R> {
32    /// Create an S-scale for loss `rho` with consistency target `0 < delta < sup ρ`.
33    ///
34    /// Rejects an unbounded loss (`rho_sup() == None`) with
35    /// `RobustError::UnboundedLoss`, since it cannot define a high-breakdown
36    /// scale and a `delta` outside `(0, sup ρ)` (for which the scale equation
37    /// has no root) with `RobustError::InvalidTuning`.
38    pub fn new(rho: R, delta: f64) -> Result<Self, RobustError> {
39        // Boundedness is the S-scale's defining precondition, not an
40        // optimization: with unbounded ρ the equation still has a unique root,
41        // but that root has zero breakdown, the opposite of the point.
42        let sup = rho.rho_sup().ok_or(RobustError::UnboundedLoss)?;
43        if delta.is_finite() && delta > 0.0 && delta < sup {
44            Ok(Self {
45                rho,
46                delta,
47                control: Control::default(),
48            })
49        } else {
50            Err(RobustError::InvalidTuning { value: delta })
51        }
52    }
53
54    /// Create a Fisher-consistent S-scale by setting `δ = E_Φ[ρ]`, evaluated by
55    /// Gauss–Hermite quadrature with `quad_points` nodes.
56    pub fn fisher_consistent(rho: R, quad_points: usize) -> Result<Self, RobustError> {
57        let (nodes, weights) = gauss_hermite(quad_points);
58        let delta = nodes
59            .iter()
60            .zip(&weights)
61            .map(|(&x, &w)| w * rho.rho(x))
62            .sum();
63        Self::new(rho, delta)
64    }
65
66    /// Override the default convergence control.
67    pub fn with_control(mut self, control: Control) -> Self {
68        self.control = control;
69        self
70    }
71
72    /// The consistency target `δ`.
73    pub fn delta(&self) -> f64 {
74        self.delta
75    }
76}
77
78impl<R: RhoFunction> ScaleEstimator for SScale<R> {
79    fn scale(&self, residuals: &[f64]) -> Result<Scale, RobustError> {
80        let n = residuals.len();
81        if n == 0 {
82            return Err(RobustError::InsufficientData { needed: 1, got: 0 });
83        }
84
85        // MAD start (about the median); any positive scale seeds the fixed point.
86        let mut buf = residuals.to_vec();
87        let med = median(&mut buf);
88        for (b, &r) in buf.iter_mut().zip(residuals) {
89            *b = (r - med).abs();
90        }
91        let mut s = 1.482_602_218_505_602 * median(&mut buf);
92        if !(s.is_finite() && s > 0.0) {
93            return Err(RobustError::DegenerateScale); // ≥ half the residuals tied
94        }
95
96        // Fixed point on (1/n) Σ ρ(rᵢ/s) = δ:  s ← s·√(mean ρ(rᵢ/s) / δ).
97        for _ in 0..self.control.max_iter {
98            let mean_rho = residuals.iter().map(|&r| self.rho.rho(r / s)).sum::<f64>() / n as f64;
99            let s_next = s * (mean_rho / self.delta).sqrt();
100            if !(s_next.is_finite() && s_next > 0.0) {
101                return Err(RobustError::DegenerateScale);
102            }
103            if (s_next - s).abs() <= self.control.tol * s {
104                return Scale::new(s_next);
105            }
106            s = s_next;
107        }
108        Err(RobustError::NonConvergence {
109            iters: self.control.max_iter,
110        })
111    }
112}
113
114/// Median via in-place total-order sort.
115fn median(v: &mut [f64]) -> f64 {
116    v.sort_unstable_by(f64::total_cmp);
117    let n = v.len();
118    let mid = n / 2;
119    if n % 2 == 1 {
120        v[mid]
121    } else {
122        0.5 * (v[mid - 1] + v[mid])
123    }
124}
robust_rs_core/scale/s_scale.rs

robust_rs_core/scale/
s_scale.rs
