robust_rs/regression/
theil_sen.rs

1//! The Theil–Sen estimator (Theil 1950; Sen 1968): the median of pairwise
2//! slopes.
3//!
4//! For simple (single-predictor) linear regression `y ≈ a + b·x`, the slope is
5//! the median of `(yⱼ − yᵢ)/(xⱼ − xᵢ)` over all pairs `i < j` and the intercept
6//! is `median(yᵢ − b·xᵢ)`. Pairs with `xᵢ = xⱼ` give an undefined slope and are
7//! **skipped**, so the slope is the median over the `n(n−1)/2` pairs *minus the
8//! x-ties*. High breakdown (`1 − 1/√2 ≈ 0.293`) with no tuning constant.
9//!
10//! Like [`crate::regression::LtsFit`] this is a bespoke result type, not an
11//! M-estimator. Theil–Sen has a known Gaussian efficiency, but it is not
12//! `ρ`-derivable, so it is not reported through the `ρ`-based
13//! [`crate::estimator::RobustEstimator`] API; only the breakdown point (a settled
14//! constant) is exposed.
15
16use crate::util::median;
17use robust_rs_core::error::RobustError;
18
19/// Theil–Sen asymptotic breakdown point, `1 − 1/√2 ≈ 0.293`.
20pub const THEIL_SEN_BREAKDOWN: f64 = 1.0 - std::f64::consts::FRAC_1_SQRT_2;
21
22/// A fitted Theil–Sen simple-regression line `y ≈ intercept + slope·x`.
23#[derive(Debug, Clone, Copy, PartialEq)]
24pub struct TheilSenFit {
25    /// The median-of-pairwise-slopes estimate.
26    pub slope: f64,
27    /// The intercept, `median(yᵢ − slope·xᵢ)`.
28    pub intercept: f64,
29}
30
31impl TheilSenFit {
32    /// The estimated slope.
33    pub fn slope(&self) -> f64 {
34        self.slope
35    }
36    /// The estimated intercept.
37    pub fn intercept(&self) -> f64 {
38        self.intercept
39    }
40    /// The predicted response at `x`.
41    pub fn predict(&self, x: f64) -> f64 {
42        self.intercept + self.slope * x
43    }
44    /// The asymptotic breakdown point (`THEIL_SEN_BREAKDOWN` = `1 − 1/√2`).
45    pub fn breakdown_point(&self) -> f64 {
46        THEIL_SEN_BREAKDOWN
47    }
48}
49
50/// Fit a Theil–Sen line to paired `(x, y)` data.
51///
52/// Returns [`RobustError::DimensionMismatch`] if the lengths differ,
53/// [`RobustError::InsufficientData`] for fewer than two points and
54/// [`RobustError::SingularDesign`] if every `xᵢ` is equal (the slope is then
55/// unidentifiable; every pair is an x-tie).
56pub fn theil_sen(x: &[f64], y: &[f64]) -> Result<TheilSenFit, RobustError> {
57    let n = x.len();
58    if y.len() != n {
59        return Err(RobustError::DimensionMismatch {
60            expected: n,
61            got: y.len(),
62        });
63    }
64    if n < 2 {
65        return Err(RobustError::InsufficientData { needed: 2, got: n });
66    }
67
68    // Pairwise slopes, skipping x-ties (an undefined slope). Comparing `dx` to
69    // the exact literal 0.0 is the intended test, a true tie, not a tolerance.
70    let mut slopes = Vec::with_capacity(n * (n - 1) / 2);
71    for i in 0..n {
72        for j in (i + 1)..n {
73            let dx = x[j] - x[i];
74            if dx != 0.0 {
75                slopes.push((y[j] - y[i]) / dx);
76            }
77        }
78    }
79    if slopes.is_empty() {
80        return Err(RobustError::SingularDesign); // all x tied
81    }
82    let slope = median(&mut slopes);
83
84    let mut offsets: Vec<f64> = x.iter().zip(y).map(|(&xi, &yi)| yi - slope * xi).collect();
85    let intercept = median(&mut offsets);
86
87    Ok(TheilSenFit { slope, intercept })
88}
robust_rs/regression/theil_sen.rs

robust_rs/regression/
theil_sen.rs
