solow_stats/
srange.rs

1//! The studentized range distribution used by Tukey's HSD.
2//!
3//! The studentized range `Q = W / U`, where `W` is the range of `k` i.i.d.
4//! standard-normal variates and `νU² ~ χ²_ν` is an independent scaled
5//! chi-squared, has CDF
6//!
7//! ```text
8//! P(Q ≤ q) = ∫₀^∞ f_U(u) · F_W(q·u) du,
9//! F_W(w)   = k ∫_{-∞}^{∞} φ(z) [Φ(z) − Φ(z − w)]^{k−1} dz,
10//! ```
11//!
12//! where `f_U` is the density of `U = √(χ²_ν/ν)`. Both integrals are evaluated
13//! with fixed-order Gauss–Legendre quadrature; the node counts are chosen so
14//! the CDF reproduces the reference (scipy `studentized_range`) to better than
15//! `1e-12`. The quantile function inverts the CDF by bisection.
16
17use solow_distributions::special::lgamma;
18use solow_distributions::{norm_cdf, norm_pdf};
19
20/// Number of Gauss–Legendre nodes for the inner (range) integral.
21const N_INNER: usize = 100;
22/// Number of Gauss–Legendre nodes for the outer (chi) integral.
23const N_OUTER: usize = 80;
24
25/// Gauss–Legendre nodes and weights on `[-1, 1]` for `n` points.
26///
27/// Computed by Newton's method on the Legendre polynomial `P_n`, which is the
28/// standard self-contained construction (no tables).
29fn gauss_legendre(n: usize) -> (Vec<f64>, Vec<f64>) {
30    let mut x = vec![0.0; n];
31    let mut w = vec![0.0; n];
32    let m = n.div_ceil(2);
33    let nf = n as f64;
34    for i in 0..m {
35        // Initial guess for the i-th root (Chebyshev approximation).
36        let mut z = (std::f64::consts::PI * (i as f64 + 0.75) / (nf + 0.5)).cos();
37        let mut pp = 0.0;
38        for _ in 0..100 {
39            // Evaluate the Legendre polynomial P_n(z) and its value P_{n-1}.
40            let mut p1 = 1.0;
41            let mut p2 = 0.0;
42            for j in 0..n {
43                let p3 = p2;
44                p2 = p1;
45                let jf = j as f64;
46                p1 = ((2.0 * jf + 1.0) * z * p2 - jf * p3) / (jf + 1.0);
47            }
48            // Derivative via the recurrence relation.
49            pp = nf * (z * p1 - p2) / (z * z - 1.0);
50            let z1 = z;
51            z = z1 - p1 / pp;
52            if (z - z1).abs() <= 1e-15 {
53                break;
54            }
55        }
56        x[i] = -z;
57        x[n - 1 - i] = z;
58        let wt = 2.0 / ((1.0 - z * z) * pp * pp);
59        w[i] = wt;
60        w[n - 1 - i] = wt;
61    }
62    (x, w)
63}
64
65/// Standard-normal CDF (re-exported for clarity).
66#[inline]
67fn phi_cdf(z: f64) -> f64 {
68    norm_cdf(z)
69}
70
71/// `F_W(w)`: CDF of the range of `k` i.i.d. standard normals at `w`.
72fn range_cdf(w: f64, k: f64, nodes: &(Vec<f64>, Vec<f64>)) -> f64 {
73    if w <= 0.0 {
74        return 0.0;
75    }
76    // The integrand `φ(z)[Φ(z)−Φ(z−w)]^{k−1}` is negligible outside this band.
77    let a = -8.0;
78    let b = 8.0 + w;
79    let half = 0.5 * (b - a);
80    let mid = 0.5 * (b + a);
81    let (x, wt) = nodes;
82    let mut acc = 0.0;
83    for i in 0..x.len() {
84        let z = half * x[i] + mid;
85        let inner = phi_cdf(z) - phi_cdf(z - w);
86        acc += wt[i] * norm_pdf(z) * inner.powf(k - 1.0);
87    }
88    k * half * acc
89}
90
91/// CDF of the studentized range with `k` groups and `df` degrees of freedom.
92pub fn srange_cdf(q: f64, k: f64, df: f64) -> f64 {
93    if q <= 0.0 {
94        return 0.0;
95    }
96    let inner = gauss_legendre(N_INNER);
97    if !df.is_finite() {
98        return range_cdf(q, k, &inner);
99    }
100    let outer = gauss_legendre(N_OUTER);
101    // log of the normalising constant of f_U(u) = c u^{ν−1} exp(−ν u²/2).
102    let logc = (df / 2.0) * df.ln() - (df / 2.0 - 1.0) * std::f64::consts::LN_2 - lgamma(df / 2.0);
103    let a = 1e-9;
104    let b = 1.0 + 10.0 / df.sqrt();
105    let half = 0.5 * (b - a);
106    let mid = 0.5 * (b + a);
107    let (x, wt) = &outer;
108    let mut acc = 0.0;
109    for i in 0..x.len() {
110        let u = half * x[i] + mid;
111        let f_u = (logc + (df - 1.0) * u.ln() - df * u * u / 2.0).exp();
112        acc += wt[i] * f_u * range_cdf(q * u, k, &inner);
113    }
114    let v = half * acc;
115    v.clamp(0.0, 1.0)
116}
117
118/// Survival function `1 − CDF` of the studentized range.
119pub fn srange_sf(q: f64, k: f64, df: f64) -> f64 {
120    1.0 - srange_cdf(q, k, df)
121}
122
123/// Quantile (inverse CDF) of the studentized range at probability `p`.
124///
125/// Found by bisection on the monotone CDF; the bracket `[0, 100]` covers every
126/// practical case.
127pub fn srange_ppf(p: f64, k: f64, df: f64) -> f64 {
128    if p <= 0.0 {
129        return 0.0;
130    }
131    if p >= 1.0 {
132        return f64::INFINITY;
133    }
134    let mut lo = 0.0;
135    let mut hi = 100.0;
136    for _ in 0..200 {
137        let mid = 0.5 * (lo + hi);
138        if srange_cdf(mid, k, df) < p {
139            lo = mid;
140        } else {
141            hi = mid;
142        }
143        if hi - lo <= 1e-12 * (1.0 + hi) {
144            break;
145        }
146    }
147    0.5 * (lo + hi)
148}
149
150#[cfg(test)]
151mod tests {
152    use super::*;
153
154    #[test]
155    fn gauss_legendre_integrates_polynomial() {
156        // ∫_{-1}^{1} x² dx = 2/3, exact for 2 nodes.
157        let (x, w) = gauss_legendre(2);
158        let val: f64 = x.iter().zip(&w).map(|(&xi, &wi)| wi * xi * xi).sum();
159        assert!((val - 2.0 / 3.0).abs() < 1e-12);
160    }
161
162    #[test]
163    fn cdf_is_monotone_and_bounded() {
164        let a = srange_cdf(2.0, 3.0, 20.0);
165        let b = srange_cdf(4.0, 3.0, 20.0);
166        assert!(a > 0.0 && a < b && b < 1.0);
167    }
168
169    #[test]
170    fn ppf_inverts_cdf() {
171        let q = srange_ppf(0.95, 4.0, 30.0);
172        let p = srange_cdf(q, 4.0, 30.0);
173        assert!((p - 0.95).abs() < 1e-8);
174    }
175}
solow_stats/srange.rs

solow_stats/
srange.rs
