solow_stats/
influence.rs

1//! Collinearity and goodness-of-fit-to-normal diagnostics.
2//!
3//! Provides the variance-inflation factor ([`variance_inflation_factor`]) for a
4//! single regressor and the Lilliefors / Kolmogorov–Smirnov test of normality
5//! ([`lilliefors`], [`kstest_normal`]).
6//!
7//! Mirrors the reference `…stats.outliers_influence.variance_inflation_factor`
8//! and `…stats.diagnostic.lilliefors` / `kstest_normal`.
9
10use ndarray::{Array1, Array2};
11use solow_core::error::{Error, Result};
12use solow_distributions::norm_cdf;
13use solow_regression::LinearModel;
14
15/// Variance-inflation factor of column `exog_idx` of the design matrix `exog`.
16///
17/// Regresses column `exog_idx` on all *other* columns (no separate constant is
18/// added — `exog` should already contain whatever constant the model uses) and
19/// returns `VIF = 1 / (1 − R²)` of that auxiliary regression. Mirrors the
20/// reference `variance_inflation_factor`.
21pub fn variance_inflation_factor(exog: &Array2<f64>, exog_idx: usize) -> Result<f64> {
22    let (n, k) = exog.dim();
23    if exog_idx >= k {
24        return Err(Error::Value("exog_idx out of range".into()));
25    }
26    if k < 2 {
27        return Err(Error::Value("need at least two columns for a VIF".into()));
28    }
29    // x_i is the target column; x_noti is everything else (column order kept).
30    let mut x_noti = Array2::<f64>::zeros((n, k - 1));
31    let mut cc = 0usize;
32    for j in 0..k {
33        if j == exog_idx {
34            continue;
35        }
36        for i in 0..n {
37            x_noti[[i, cc]] = exog[[i, j]];
38        }
39        cc += 1;
40    }
41    let x_i = exog.column(exog_idx).to_owned();
42    let res = LinearModel::ols(x_i, x_noti)?.fit()?;
43    Ok(1.0 / (1.0 - res.rsquared))
44}
45
46/// Assumed reference distribution for the Lilliefors test.
47#[derive(Debug, Clone, Copy, PartialEq, Eq)]
48pub enum LillieforsDist {
49    /// Normal distribution with mean and variance estimated from the sample.
50    Norm,
51}
52
53/// Two-sided Kolmogorov–Smirnov statistic of the sorted standardized sample
54/// `z` against the standard-normal CDF.
55///
56/// `D = max(D⁺, D⁻)` with `D⁺ = maxᵢ (i/n − Φ(z₍ᵢ₎))` and
57/// `D⁻ = maxᵢ (Φ(z₍ᵢ₎) − (i−1)/n)` for `i = 1 … n` (sorted ascending). This is
58/// the reference `ksstat(z, norm.cdf, 'two_sided')`.
59fn ks_stat_normal(z: &mut [f64]) -> f64 {
60    z.sort_by(|a, b| a.total_cmp(b));
61    let n = z.len() as f64;
62    let mut d_plus = f64::NEG_INFINITY;
63    let mut d_min = f64::NEG_INFINITY;
64    for (idx, &zi) in z.iter().enumerate() {
65        let cdf = norm_cdf(zi);
66        let i = idx as f64; // 0-based
67        let dp = (i + 1.0) / n - cdf; // (i+1)/n - F
68        let dm = cdf - i / n; // F - i/n
69        if dp > d_plus {
70            d_plus = dp;
71        }
72        if dm > d_min {
73            d_min = dm;
74        }
75    }
76    d_plus.max(d_min)
77}
78
79/// Dalal–Wilkinson approximation of the Lilliefors p-value for normality.
80///
81/// Valid (per the reference) for p-values below ~0.1; this is the closed-form
82/// `pval_lf` used by the reference's `pvalmethod="approx"`. For `n > 100` the
83/// statistic is rescaled and `n` capped at 100.
84fn pval_lf(d_max: f64, n: usize) -> f64 {
85    let mut d = d_max;
86    let mut nn = n as f64;
87    if n > 100 {
88        d *= (nn / 100.0).powf(0.49);
89        nn = 100.0;
90    }
91    (-7.01256 * d * d * (nn + 2.78019) + 2.99587 * d * (nn + 2.78019).sqrt() - 0.122119
92        + 0.974598 / nn.sqrt()
93        + 1.67997 / nn)
94        .exp()
95}
96
97/// Lilliefors test of normality with estimated parameters.
98///
99/// Standardizes `x` by its sample mean and (ddof = 1) standard deviation, then
100/// returns the two-sided KS statistic against the standard normal together with
101/// the Dalal–Wilkinson approximate p-value (the reference's
102/// `pvalmethod="approx"`). The statistic is closed form; the p-value is the
103/// closed-form approximation. Requires `n ≥ 4`.
104///
105/// Mirrors the reference `lilliefors(x, dist="norm", pvalmethod="approx")`.
106pub fn lilliefors(x: &Array1<f64>, dist: LillieforsDist) -> Result<(f64, f64)> {
107    let LillieforsDist::Norm = dist;
108    let n = x.len();
109    if n < 4 {
110        return Err(Error::Value(
111            "Lilliefors test requires at least 4 observations".into(),
112        ));
113    }
114    let nf = n as f64;
115    let mean = x.sum() / nf;
116    // Sample standard deviation (ddof = 1).
117    let var = x.iter().map(|&v| (v - mean) * (v - mean)).sum::<f64>() / (nf - 1.0);
118    let sd = var.sqrt();
119    let mut z: Vec<f64> = x.iter().map(|&v| (v - mean) / sd).collect();
120    let d_ks = ks_stat_normal(&mut z);
121    let pval = pval_lf(d_ks, n);
122    Ok((d_ks, pval))
123}
124
125/// Alias for the Lilliefors normality test, matching the reference name
126/// `kstest_normal`. Returns `(statistic, pvalue)`.
127pub fn kstest_normal(x: &Array1<f64>) -> Result<(f64, f64)> {
128    lilliefors(x, LillieforsDist::Norm)
129}
130
131#[cfg(test)]
132mod tests {
133    use super::*;
134    use ndarray::array;
135
136    #[test]
137    fn vif_orthogonal_columns_is_one() {
138        // Two orthogonal slopes plus a constant: VIF of a slope == 1.
139        let x = array![
140            [1.0, 1.0, 1.0],
141            [1.0, 1.0, -1.0],
142            [1.0, -1.0, 1.0],
143            [1.0, -1.0, -1.0],
144        ];
145        let v = variance_inflation_factor(&x, 1).unwrap();
146        assert!((v - 1.0).abs() < 1e-9);
147    }
148
149    #[test]
150    fn lilliefors_statistic_in_unit_range() {
151        let x = array![0.1, -0.5, 0.3, 1.2, -0.7, 0.4, -0.2, 0.9, -1.1, 0.05];
152        let (stat, p) = lilliefors(&x, LillieforsDist::Norm).unwrap();
153        assert!((0.0..=1.0).contains(&stat));
154        assert!(p > 0.0);
155        let (stat2, _) = kstest_normal(&x).unwrap();
156        assert!((stat - stat2).abs() < 1e-15);
157    }
158}
solow_stats/influence.rs

solow_stats/
influence.rs
