robust_rs_core/scale/
huber_proposal2.rs

1//! Huber's "Proposal 2" simultaneous location–scale estimate.
2
3use crate::error::RobustError;
4use crate::scale::ScaleEstimator;
5use crate::types::Scale;
6
7/// Huber's Proposal 2 scale: solve `(1/n) Σ ψ((rᵢ−μ)/s)² = β` (with
8/// `β = E_Φ[ψ²]`) jointly with the location, by fixed-point iteration.
9#[derive(Debug, Clone, Copy)]
10pub struct HuberProposal2 {
11    /// Tuning constant `k` for the ψ used in the scale equation.
12    pub k: f64,
13}
14
15impl Default for HuberProposal2 {
16    fn default() -> Self {
17        Self { k: 1.345 }
18    }
19}
20
21impl ScaleEstimator for HuberProposal2 {
22    fn scale(&self, residuals: &[f64]) -> Result<Scale, RobustError> {
23        let n = residuals.len();
24        if n == 0 {
25            return Err(RobustError::InsufficientData { needed: 1, got: 0 });
26        }
27        let k = self.k;
28        if !(k.is_finite() && k > 0.0) {
29            return Err(RobustError::InvalidTuning { value: k }); // field is public/unvalidated
30        }
31
32        // Consistency constant β = E_Φ[ψ_k(Z)²], in CLOSED FORM. Do not compute
33        // this by Gauss–Hermite quadrature: ψ² has a kink at ±k, so GH converges
34        // at ~O(1/n) and oscillates (~0.3% error even at 128 nodes), which would
35        // bias every scale it produces. The closed form is exact.
36        let phi_k = (-0.5 * k * k).exp() / (2.0 * std::f64::consts::PI).sqrt();
37        let cdf_k = normal_cdf(k);
38        let beta = (2.0 * cdf_k - 1.0 - 2.0 * k * phi_k) + 2.0 * k * k * (1.0 - cdf_k);
39
40        // Robust start: median and MAD.
41        let mut buf = residuals.to_vec();
42        let mut mu = median(&mut buf);
43        for (b, &r) in buf.iter_mut().zip(residuals) {
44            *b = (r - mu).abs();
45        }
46        let mut s = 1.482_602_218_505_602 * median(&mut buf);
47        if !(s.is_finite() && s > 0.0) {
48            return Err(RobustError::DegenerateScale); // ≥ half the residuals tied
49        }
50
51        // No `Control` is threaded through the ScaleEstimator trait, so mirror
52        // the crate's default stopping rule.
53        const TOL: f64 = 1e-8;
54        const MAX_ITER: usize = 100;
55
56        // Fixed point on the joint location–scale equations:
57        //   winsorize each rᵢ to [μ−ks, μ+ks];  μ ← mean(winsorized);
58        //   s ← sqrt( Σ (winsorized − μ)² / (n·β) ).
59        for _ in 0..MAX_ITER {
60            let (lo, hi) = (mu - k * s, mu + k * s);
61
62            let mut sum = 0.0_f64;
63            for &r in residuals {
64                sum += r.clamp(lo, hi);
65            }
66            let mu1 = sum / n as f64;
67
68            let mut ss = 0.0_f64;
69            for &r in residuals {
70                let d = r.clamp(lo, hi) - mu1;
71                ss += d * d;
72            }
73            let s1 = (ss / (n as f64 * beta)).sqrt();
74
75            if !(s1.is_finite() && s1 > 0.0) {
76                return Err(RobustError::DegenerateScale);
77            }
78            if (mu - mu1).abs() < TOL * s && (s - s1).abs() < TOL * s {
79                return Scale::new(s1);
80            }
81            mu = mu1;
82            s = s1;
83        }
84        Err(RobustError::NonConvergence { iters: MAX_ITER })
85    }
86}
87
88/// Median via in-place total-order sort.
89fn median(v: &mut [f64]) -> f64 {
90    v.sort_unstable_by(f64::total_cmp);
91    let n = v.len();
92    let mid = n / 2;
93    if n % 2 == 1 {
94        v[mid]
95    } else {
96        0.5 * (v[mid - 1] + v[mid])
97    }
98}
99
100/// Standard normal CDF via `libm::erf` (full double precision), keeping the core
101/// crate's special-function surface dependency-light without `statrs`.
102fn normal_cdf(x: f64) -> f64 {
103    0.5 * (1.0 + libm::erf(x / std::f64::consts::SQRT_2))
104}
robust_rs_core/scale/huber_proposal2.rs

robust_rs_core/scale/
huber_proposal2.rs
