regression-diagnostics 0.2.0

use ndarray::Array1;

use crate::OlsFit;

/// Internally studentized residuals `rᵢ = eᵢ / (s·√(1 − hᵢ))`.
///
/// Each raw residual is scaled by its own standard deviation, which shrinks
/// toward zero as leverage `hᵢ` grows. "Internally" means the global residual
/// standard error `s` (which was itself computed *using* observation `i`) is
/// used. This is R's `rstandard`.
pub fn internally_studentized_residuals(fit: &OlsFit) -> Array1<f64> {
    let s = fit.residual_standard_error();
    let lev = fit.leverage();
    Array1::from_shape_fn(fit.n_observations(), |i| {
        let denom = s * (1.0 - lev[i]).max(0.0).sqrt();
        if denom > 0.0 {
            fit.residuals()[i] / denom
        } else {
            f64::NAN
        }
    })
}

/// Externally studentized residuals `tᵢ = eᵢ / (s₍ᵢ₎·√(1 − hᵢ))`.
///
/// "Externally" means observation `i` is *excluded* from the variance estimate
/// used to scale its own residual:
///
/// `s₍ᵢ₎² = [(n − p)·s² − eᵢ²/(1 − hᵢ)] / (n − p − 1)`
///
/// so a genuine outlier no longer inflates the very scale it is judged against —
/// which is why this form (R's `rstudent`) is the one that follows a Student-t
/// distribution and underlies DFFITS. Requires `n − p − 1 ≥ 1`; entries are
/// `NaN` otherwise.
pub fn externally_studentized_residuals(fit: &OlsFit) -> Array1<f64> {
    let n = fit.n_observations() as f64;
    let p = fit.n_parameters() as f64;
    let df = n - p; // residual df of the full fit
    let s2 = fit.residual_variance();
    let lev = fit.leverage();

    Array1::from_shape_fn(fit.n_observations(), |i| {
        let h = lev[i];
        let one_minus_h = (1.0 - h).max(0.0);
        if df - 1.0 < 1.0 || one_minus_h <= 0.0 {
            return f64::NAN;
        }
        let e = fit.residuals()[i];
        let s2_i = (df * s2 - e * e / one_minus_h) / (df - 1.0);
        if s2_i <= 0.0 {
            return f64::NAN;
        }
        e / (s2_i.sqrt() * one_minus_h.sqrt())
    })
}