qualia_core_db/solvers/statistics/
regression.rs

1//! Regression — ordinary least squares. Simple (one-predictor) linear regression
2//! with the full inferential output: coefficients, R², residual standard error, and
3//! real t-based standard errors / p-values from the [`distributions`](super::distributions)
4//! library (no placeholder significance).
5//!
6//! Multiple linear regression (the normal-equations / QR solve over
7//! `solvers::linear_algebra`) is the natural next module here; simple OLS is a
8//! complete capability on its own and is what the domain libs need first.
9
10use super::descriptive::{covariance, mean, variance};
11use super::distributions::students_t;
12
13/// Ordinary-least-squares fit of `y = intercept + slope·x`.
14#[derive(Debug, Clone, Copy, PartialEq)]
15pub struct LinearRegression {
16    pub slope: f64,
17    pub intercept: f64,
18    /// Coefficient of determination R² (fraction of variance explained).
19    pub r_squared: f64,
20    /// Residual standard error `s = √(SSE/(n−2))`.
21    pub residual_std_error: f64,
22    pub slope_std_error: f64,
23    pub slope_t: f64,
24    /// Two-sided p-value for `slope = 0` (df = n−2).
25    pub slope_p_value: f64,
26    pub intercept_std_error: f64,
27    pub intercept_p_value: f64,
28    pub n: usize,
29}
30
31/// Simple linear regression of `y` on `x`. `None` if the lengths differ, `n < 3`
32/// (need `n−2 ≥ 1` residual degrees of freedom for inference), or `x` has zero
33/// variance (slope undefined).
34pub fn simple_linear_regression(x: &[f64], y: &[f64]) -> Option<LinearRegression> {
35    let n = x.len();
36    if n != y.len() || n < 3 {
37        return None;
38    }
39    let mx = mean(x)?;
40    let my = mean(y)?;
41    let var_x = variance(x, false)?; // population moment is fine; ratios cancel n
42    if var_x <= 0.0 {
43        return None;
44    }
45    let cov_xy = covariance(x, y, false)?;
46    let slope = cov_xy / var_x;
47    let intercept = my - slope * mx;
48
49    // Sums of squares.
50    let mut ss_tot = 0.0; // Σ(y-ȳ)²
51    let mut ss_res = 0.0; // Σ(y-ŷ)²
52    let mut sxx = 0.0; // Σ(x-x̄)²
53    for i in 0..n {
54        let yhat = intercept + slope * x[i];
55        ss_tot += (y[i] - my).powi(2);
56        ss_res += (y[i] - yhat).powi(2);
57        sxx += (x[i] - mx).powi(2);
58    }
59
60    let df = (n - 2) as f64;
61    let r_squared = if ss_tot > 0.0 {
62        1.0 - ss_res / ss_tot
63    } else {
64        1.0
65    };
66    let s2 = ss_res / df; // residual variance
67    let residual_std_error = s2.sqrt();
68
69    let slope_std_error = (s2 / sxx).sqrt();
70    let (slope_t, slope_p_value) = if slope_std_error > 0.0 {
71        let t = slope / slope_std_error;
72        (t, students_t::two_sided_p(t, df))
73    } else {
74        // Perfect fit (zero residuals): slope is exact → infinitely significant.
75        (f64::INFINITY.copysign(slope), 0.0)
76    };
77
78    let intercept_std_error = (s2 * (1.0 / n as f64 + mx * mx / sxx)).sqrt();
79    let intercept_p_value = if intercept_std_error > 0.0 {
80        students_t::two_sided_p(intercept / intercept_std_error, df)
81    } else {
82        0.0
83    };
84
85    Some(LinearRegression {
86        slope,
87        intercept,
88        r_squared,
89        residual_std_error,
90        slope_std_error,
91        slope_t,
92        slope_p_value,
93        intercept_std_error,
94        intercept_p_value,
95        n,
96    })
97}
98
99#[cfg(test)]
100mod tests {
101    use super::*;
102
103    #[test]
104    fn exact_line_is_recovered() {
105        // y = 3 + 2x exactly.
106        let x = [0.0, 1.0, 2.0, 3.0, 4.0];
107        let y = [3.0, 5.0, 7.0, 9.0, 11.0];
108        let r = simple_linear_regression(&x, &y).unwrap();
109        assert!((r.slope - 2.0).abs() < 1e-9);
110        assert!((r.intercept - 3.0).abs() < 1e-9);
111        assert!((r.r_squared - 1.0).abs() < 1e-12);
112        assert!(r.residual_std_error < 1e-9);
113        assert_eq!(r.slope_p_value, 0.0); // perfect fit → exact
114    }
115
116    #[test]
117    fn noisy_trend_matches_known_fit() {
118        // Classic small dataset; OLS slope/intercept are standard reference values.
119        let x = [1.0, 2.0, 3.0, 4.0, 5.0];
120        let y = [2.1, 3.9, 6.1, 7.9, 10.2];
121        let r = simple_linear_regression(&x, &y).unwrap();
122        // Hand/np.polyfit: slope ≈ 2.0, intercept ≈ 0.02.
123        assert!((r.slope - 2.0).abs() < 0.05, "slope={}", r.slope);
124        assert!(r.intercept.abs() < 0.2, "intercept={}", r.intercept);
125        assert!(r.r_squared > 0.99, "r2={}", r.r_squared);
126        assert!(r.slope_p_value < 1e-4, "strong trend must be significant");
127    }
128
129    #[test]
130    fn no_relationship_is_not_significant() {
131        let x = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
132        let y = [5.0, 4.0, 6.0, 5.0, 6.0, 4.0]; // flat-ish noise
133        let r = simple_linear_regression(&x, &y).unwrap();
134        assert!(
135            r.slope_p_value > 0.2,
136            "no real trend: p={}",
137            r.slope_p_value
138        );
139        assert!(r.r_squared < 0.3);
140    }
141
142    #[test]
143    fn guards_degenerate_input() {
144        assert!(simple_linear_regression(&[1.0, 2.0], &[1.0, 2.0]).is_none()); // n<3
145        assert!(simple_linear_regression(&[2.0, 2.0, 2.0], &[1.0, 2.0, 3.0]).is_none()); // zero var x
146        assert!(simple_linear_regression(&[1.0, 2.0, 3.0], &[1.0, 2.0]).is_none());
147        // mismatch
148    }
149}
qualia_core_db/solvers/statistics/regression.rs

qualia_core_db/solvers/statistics/
regression.rs
