solow-stats 0.7.3

//! Noncentral-t distribution CDF, used by one-sample t-test power.
//!
//! Implements Lenth's algorithm (Applied Statistics AS 243, 1989) for the CDF
//! of the noncentral t-distribution with `df` degrees of freedom and
//! noncentrality `nc`. The series reproduces the reference (`scipy.special`'s
//! `nctdtr`) to better than `1e-12` across the parameter ranges relevant to
//! power calculations.

use solow_distributions::norm_cdf;
use solow_distributions::special::{betainc, lgamma};

const TWO_OVER_PI_SQRT: f64 = 0.797_884_560_802_865_4; // sqrt(2/pi)

/// CDF of the noncentral t-distribution: `P(T ≤ t)` with `df` d.o.f. and
/// noncentrality `nc`.
pub fn nct_cdf(t: f64, df: f64, nc: f64) -> f64 {
    // Use the reflection identity F(t; df, nc) = 1 − F(−t; df, −nc) so the
    // series is always evaluated for a non-negative argument.
    let (tt, dl, negdel) = if t < 0.0 {
        (-t, -nc, true)
    } else {
        (t, nc, false)
    };

    let x = tt * tt / (tt * tt + df);
    if x <= 0.0 {
        let res = norm_cdf(-dl);
        return if negdel { 1.0 - res } else { res };
    }

    let lambda = dl * dl;
    let mut p = 0.5 * (-0.5 * lambda).exp();
    let mut q = TWO_OVER_PI_SQRT * p * dl;
    let mut s = 0.5 - p;
    let mut a = 0.5;
    let b = 0.5 * df;
    let rxb = (1.0 - x).powf(b);
    let albeta = lgamma(a) + lgamma(b) - lgamma(a + b);

    let mut xodd = betainc(a, b, x);
    let mut godd = 2.0 * rxb * (a * x.ln() - albeta).exp();
    let mut xeven = 1.0 - rxb;
    let mut geven = b * x * rxb;
    let mut tnc = p * xodd + q * xeven;

    let errbd = 1e-14;
    let maxit = 2000;
    let mut it = 1;
    while it <= maxit {
        a += 1.0;
        xodd -= godd;
        xeven -= geven;
        godd *= x * (a + b - 1.0) / a;
        geven *= x * (a + b - 0.5) / (a + 0.5);
        p *= lambda / (2.0 * it as f64);
        q *= lambda / (2.0 * it as f64 + 1.0);
        s -= p;
        tnc += p * xodd + q * xeven;
        let err = 2.0 * s * (xodd - godd);
        if err.abs() < errbd {
            break;
        }
        it += 1;
    }

    tnc += norm_cdf(-dl);
    let tnc = tnc.clamp(0.0, 1.0);
    if negdel {
        1.0 - tnc
    } else {
        tnc
    }
}

/// Survival function `1 − CDF` of the noncentral t-distribution.
pub fn nct_sf(t: f64, df: f64, nc: f64) -> f64 {
    1.0 - nct_cdf(t, df, nc)
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn central_reduces_to_t() {
        // With nc = 0 the noncentral t reduces to Student's t.
        let got = nct_cdf(1.5, 10.0, 0.0);
        let want = solow_distributions::t_cdf(1.5, 10.0);
        assert!((got - want).abs() < 1e-10, "{got} vs {want}");
    }

    #[test]
    fn reflection_identity() {
        let a = nct_cdf(2.0, 12.0, 1.0);
        let b = 1.0 - nct_cdf(-2.0, 12.0, -1.0);
        assert!((a - b).abs() < 1e-12);
    }
}