ggplot_rs/stat/
dist.rs

1//! Student-t quantile used for confidence-interval multipliers.
2//!
3//! Per this crate's rule (no statistical calculation lives here), the quantile
4//! itself is computed by anofox-regression (`statrs::StudentsT`) when the
5//! `regression` feature is enabled. Without that feature there is no
6//! distribution library available, so we fall back to the large-sample normal
7//! 97.5% quantile `1.96` — enable `regression` for exact `qt(p, df)`.
8
9/// Standard-normal 0.975 quantile — the large-sample CI multiplier used as the
10/// fallback when the `regression` feature (and thus a t-distribution) is absent.
11pub(crate) const Z_975: f64 = 1.959_963_984_540_054;
12
13/// Student-t quantile: the `t` with `P(T ≤ t) = p` for `df` degrees of freedom
14/// (R's `qt(p, df)`). `df` may be non-integer (e.g. loess effective df).
15///
16/// With the `regression` feature this delegates to anofox-regression's
17/// `StudentsT` (statrs) for an exact value; without it, it returns the normal
18/// approximation [`Z_975`] (correct only for `p = 0.975`).
19#[cfg(feature = "regression")]
20pub fn qt(p: f64, df: f64) -> f64 {
21    use anofox_regression::distributions::{ContinuousCDF, StudentsT};
22    if df <= 0.0 {
23        return f64::NAN;
24    }
25    StudentsT::new(0.0, 1.0, df)
26        .map(|d| d.inverse_cdf(p))
27        .unwrap_or(Z_975)
28}
29
30/// See the `regression`-enabled variant above; without a distribution library
31/// this returns the large-sample normal 0.975 quantile.
32#[cfg(not(feature = "regression"))]
33pub fn qt(_p: f64, _df: f64) -> f64 {
34    Z_975
35}
36
37/// Standard-normal quantile (probit), R's `qnorm(p)`.
38///
39/// With the `regression` feature this delegates to anofox-regression's `Normal`
40/// (statrs) for an exact value; without it, it uses the Abramowitz & Stegun
41/// 26.2.23 rational approximation (|error| < ~4.5e-4 in z).
42#[cfg(feature = "regression")]
43pub fn qnorm(p: f64) -> f64 {
44    use anofox_regression::distributions::{ContinuousCDF, Normal};
45    if p <= 0.0 {
46        return f64::NEG_INFINITY;
47    }
48    if p >= 1.0 {
49        return f64::INFINITY;
50    }
51    Normal::new(0.0, 1.0)
52        .map(|d| d.inverse_cdf(p))
53        .unwrap_or(0.0)
54}
55
56/// See the `regression`-enabled variant; this is the Abramowitz & Stegun
57/// rational approximation used when no distribution library is available.
58#[cfg(not(feature = "regression"))]
59pub fn qnorm(p: f64) -> f64 {
60    if p <= 0.0 {
61        return f64::NEG_INFINITY;
62    }
63    if p >= 1.0 {
64        return f64::INFINITY;
65    }
66    // A&S 26.2.23, using symmetry about p = 0.5.
67    fn rational_approx(t: f64) -> f64 {
68        let (c0, c1, c2) = (2.515_517, 0.802_853, 0.010_328);
69        let (d1, d2, d3) = (1.432_788, 0.189_269, 0.001_308);
70        t - (c0 + c1 * t + c2 * t * t) / (1.0 + d1 * t + d2 * t * t + d3 * t * t * t)
71    }
72    if p < 0.5 {
73        -rational_approx((-2.0 * p.ln()).sqrt())
74    } else if p > 0.5 {
75        rational_approx((-2.0 * (1.0 - p).ln()).sqrt())
76    } else {
77        0.0
78    }
79}
80
81/// Radius scaling for a bivariate confidence ellipse at `level` from `n` points,
82/// matching ggplot2's `stat_ellipse`: `sqrt(2 · F⁻¹(level; 2, n−1))`.
83///
84/// With the `regression` feature the F quantile comes from anofox-regression's
85/// `FisherSnedecor` (statrs). Without it, we fall back to the large-sample
86/// limit `sqrt(-2·ln(1−level))` (the χ²₂ quantile, exact as n → ∞).
87#[cfg(feature = "regression")]
88pub fn ellipse_radius(level: f64, n: usize) -> f64 {
89    use anofox_regression::distributions::{ContinuousCDF, FisherSnedecor};
90    let dfd = (n as f64 - 1.0).max(1.0);
91    match FisherSnedecor::new(2.0, dfd) {
92        Ok(f) => (2.0 * f.inverse_cdf(level)).max(0.0).sqrt(),
93        Err(_) => (-2.0 * (1.0 - level).ln()).sqrt(),
94    }
95}
96
97/// See the `regression`-enabled variant; this uses the χ²₂ closed form.
98#[cfg(not(feature = "regression"))]
99pub fn ellipse_radius(level: f64, _n: usize) -> f64 {
100    (-2.0 * (1.0 - level).ln()).sqrt()
101}
102
103#[cfg(all(test, feature = "regression"))]
104mod tests {
105    use super::{ellipse_radius, qnorm, qt};
106
107    #[test]
108    fn qnorm_matches_r() {
109        // R: qnorm(c(0.025, 0.25, 0.5, 0.75, 0.975)).
110        let cases = [
111            (0.025, -1.959_963_984_540_054),
112            (0.25, -0.674_489_750_196_082),
113            (0.5, 0.0),
114            (0.75, 0.674_489_750_196_082),
115            (0.975, 1.959_963_984_540_054),
116        ];
117        for (p, want) in cases {
118            assert!((qnorm(p) - want).abs() < 1e-9, "qnorm({p})");
119        }
120    }
121
122    #[test]
123    fn ellipse_radius_is_monotonic_and_finite() {
124        let small = ellipse_radius(0.5, 30);
125        let big = ellipse_radius(0.99, 30);
126        assert!(big > small && small > 0.0 && big.is_finite());
127    }
128
129    // Reference values from R's qt(0.975, df).
130    #[test]
131    fn qt_matches_r() {
132        let cases = [
133            (1.0, 12.706_204_736_432_1),
134            (2.0, 4.302_652_729_912_0),
135            (5.0, 2.570_581_835_636_1),
136            (10.0, 2.228_138_851_986_3),
137            (30.0, 2.042_272_456_301_4),
138            (100.0, 1.983_971_518_449_5),
139        ];
140        for (df, want) in cases {
141            let got = qt(0.975, df);
142            assert!(
143                (got - want).abs() < 1e-6,
144                "qt(0.975,{df}) = {got}, want {want}"
145            );
146        }
147    }
148}
ggplot_rs/stat/dist.rs

ggplot_rs/stat/
dist.rs
