robust_rs/location/
m_location.rs

1//! M-estimation of location by IRLS (scalar; no linear algebra needed).
2
3use robust_rs_core::error::RobustError;
4use robust_rs_core::rho::RhoFunction;
5use robust_rs_core::scale::ScaleEstimator;
6use robust_rs_core::solver::Control;
7use robust_rs_core::types::Scale;
8
9/// A fitted location estimate.
10#[derive(Debug, Clone, Copy)]
11pub struct LocationFit {
12    /// The location estimate `θ̂`.
13    pub estimate: f64,
14    /// The scale used to standardize residuals.
15    pub scale: Scale,
16    /// Iterations performed.
17    pub iters: usize,
18}
19
20/// M-estimate of location: iterate `θ ← Σ wᵢ xᵢ / Σ wᵢ` with
21/// `wᵢ = ρ.weight((xᵢ − θ)/s)` until convergence.
22pub fn m_location(
23    data: &[f64],
24    rho: &dyn RhoFunction,
25    scale: &dyn ScaleEstimator,
26    ctrl: &Control,
27) -> Result<LocationFit, RobustError> {
28    if data.is_empty() {
29        return Err(RobustError::InsufficientData { needed: 1, got: 0 });
30    }
31
32    // Scale computed once (about the data's median) and held fixed.
33    // scale.scale already returns DegenerateScale if it collapses to zero.
34    let scale_est = scale.scale(data)?;
35    let s: f64 = scale_est.get();
36
37    let mut buf = data.to_vec();
38    let mut theta = median(&mut buf); // init at the median
39
40    for i in 1..=ctrl.max_iter {
41        // wᵢ = ρ.weight((xᵢ − θ)/s);  θ ← Σ wᵢ xᵢ / Σ wᵢ
42        let (mut sw, mut swx) = (0.0_f64, 0.0_f64);
43        for &x in data {
44            let w = rho.weight((x - theta) / s);
45            sw += w;
46            swx += w * x;
47        }
48        let next = swx / sw;
49        if !next.is_finite() {
50            // 0/0 (a redescending ρ rejected every point) or non-finite data.
51            return Err(RobustError::SingularDesign);
52        }
53        // tol read as relative-to-scale: θ can legitimately be ≈0, which makes a
54        // relative-to-θ test degenerate, whereas s is a natural non-zero yardstick.
55        if (next - theta).abs() <= ctrl.tol * s {
56            return Ok(LocationFit {
57                estimate: next,
58                scale: scale_est,
59                iters: i,
60            });
61        }
62        theta = next;
63    }
64    Err(RobustError::NonConvergence {
65        iters: ctrl.max_iter,
66    })
67}
68
69/// Median via in-place total-order sort.
70fn median(v: &mut [f64]) -> f64 {
71    v.sort_unstable_by(f64::total_cmp);
72    let n = v.len();
73    let mid = n / 2;
74    if n % 2 == 1 {
75        v[mid]
76    } else {
77        0.5 * (v[mid - 1] + v[mid])
78    }
79}
robust_rs/location/m_location.rs

robust_rs/location/
m_location.rs
