robust_rs/
multivariate.rs

1//! Multivariate robust statistics: robust estimates of multivariate location
2//! and scatter and the outlier detection built on them.
3//!
4//! The organizing object here is a **location–scatter pair** `(μ̂, Σ̂)` (a
5//! robust centre and a robust covariance) from which a robust Mahalanobis
6//! distance `dᵢ = √((xᵢ − μ̂)ᵀ Σ̂⁻¹ (xᵢ − μ̂))` flags multivariate outliers, just
7//! as robust residuals flag them in regression. Four estimators produce such a
8//! pair, trading efficiency against breakdown and equivariance:
9//!
10//! - [`Mcd`]: the Minimum Covariance Determinant (Rousseeuw 1985) via FAST-MCD
11//!   (Rousseeuw & Van Driessen 1999): 50%-breakdown, affine-equivariant, the
12//!   multivariate analogue of LTS.
13//! - [`Ogk`]: the Orthogonalized Gnanadesikan–Kettenring estimator
14//!   (Maronna & Zamar 2002): a fast, deterministic, positive-definite pairwise
15//!   estimator (orthogonally, not fully affine, equivariant).
16//! - [`MScatter`]: a monotone M-estimator of location and scatter
17//!   (Maronna 1976): the direct multivariate analogue of the regression
18//!   M-estimator, reusing a [`robust_rs_core::rho::RhoFunction`] weight.
19//! - [`Tyler`]: Tyler's (1987) distribution-free M-estimator of *shape*,
20//!   normalized to unit determinant.
21//!
22//! [`mahalanobis`] exposes the distance/outlier map over *any* `(μ̂, Σ̂)` pair,
23//! together with the classical (non-robust) mean/covariance baseline.
24
25mod chi2;
26mod correction;
27mod linalg;
28
29mod m_scatter;
30pub mod mahalanobis;
31mod mcd;
32mod ogk;
33mod tyler;
34
35pub use self::m_scatter::MScatter;
36pub use self::mcd::{Mcd, McdFit};
37pub use self::ogk::Ogk;
38pub use self::tyler::{Tyler, TylerFit};
39
40use ndarray::{Array1, Array2};
41use robust_rs_core::error::RobustError;
42
43use self::chi2::chi2_quantile;
44use self::linalg::{mahalanobis_sq, spd_inverse_logdet};
45
46/// A fitted robust location–scatter estimate.
47///
48/// The shared result type of [`MScatter`] and [`Ogk`]. FAST-MCD carries extra
49/// structure (a raw estimate and the retained subset), so it returns its own
50/// [`McdFit`], but that type implements [`RobustScatter`] too, so all three are
51/// interchangeable behind the trait and the Mahalanobis/outlier map is written
52/// once against `(μ̂, Σ̂)`. [`Tyler`] does *not* use this type: it returns a
53/// bespoke [`TylerFit`] implementing no shared trait, because its shape-only
54/// distances are not `χ²`-calibrated, so this type's Gaussian outlier map would
55/// misapply to them.
56#[derive(Debug, Clone)]
57pub struct ScatterFit {
58    /// Robust location estimate `μ̂` (a `p`-vector).
59    pub location: Array1<f64>,
60    /// Robust scatter (covariance) estimate `Σ̂` (`p × p`, symmetric positive
61    /// definite).
62    pub scatter: Array2<f64>,
63    /// Robust Mahalanobis distances `dᵢ` of the rows the fit was computed on.
64    pub distances: Array1<f64>,
65    /// Final per-observation weights (their meaning is estimator-specific: a
66    /// hard `0/1` acceptance for a reweighting estimator, a smooth `w(dᵢ)` for
67    /// an M-estimator).
68    pub weights: Array1<f64>,
69}
70
71/// Quantities every fitted robust covariance estimator can report, the
72/// multivariate counterpart of [`crate::estimator::RobustEstimator`].
73pub trait RobustScatter {
74    /// Robust location estimate `μ̂`.
75    fn location(&self) -> &Array1<f64>;
76    /// Robust scatter estimate `Σ̂`.
77    fn scatter(&self) -> &Array2<f64>;
78    /// Robust Mahalanobis distances of the fitted rows.
79    fn distances(&self) -> &Array1<f64>;
80
81    /// The distance cutoff `√χ²_{p, quantile}` (e.g. `quantile = 0.975`) against
82    /// which robust Mahalanobis distances are compared to flag outliers.
83    fn distance_cutoff(&self, quantile: f64) -> f64 {
84        let p = self.scatter().nrows() as f64;
85        chi2_quantile(quantile, p).sqrt()
86    }
87
88    /// Flag each fitted observation as an outlier when its robust distance
89    /// exceeds `distance_cutoff(quantile)`. Under a `p`-variate Gaussian the
90    /// squared distances are `χ²_p`, so `quantile = 0.975` flags ≈ 2.5% of clean
91    /// data.
92    fn outliers(&self, quantile: f64) -> Vec<bool> {
93        let cut = self.distance_cutoff(quantile);
94        self.distances().iter().map(|&d| d > cut).collect()
95    }
96}
97
98impl RobustScatter for ScatterFit {
99    fn location(&self) -> &Array1<f64> {
100        &self.location
101    }
102    fn scatter(&self) -> &Array2<f64> {
103        &self.scatter
104    }
105    fn distances(&self) -> &Array1<f64> {
106        &self.distances
107    }
108}
109
110/// Robust Mahalanobis distances `dᵢ = √((xᵢ − loc)ᵀ scatter⁻¹ (xᵢ − loc))` for
111/// every row of `x`. Shared by the estimators to populate [`ScatterFit`] and
112/// re-exported publicly as [`mahalanobis::mahalanobis_distances`]. Errors if the
113/// scatter is not positive definite.
114pub(crate) fn distances_from(
115    x: &Array2<f64>,
116    loc: &Array1<f64>,
117    scatter: &Array2<f64>,
118) -> Result<Array1<f64>, RobustError> {
119    let (inv, _logdet) = spd_inverse_logdet(scatter)?;
120    Ok(mahalanobis_sq(x, loc, &inv).mapv(f64::sqrt))
121}
robust_rs/multivariate.rs

robust_rs/
multivariate.rs
