stats_claw/distributions/sampling/
students_t.rs

1//! Student's t-distribution numerics, for the
2//! [`TDistribution`].
3//!
4//! Equivalent to `scipy.stats.t(degrees_of_freedom)`. The density uses `ln_gamma`
5//! for stability, the CDF is expressed through the regularized incomplete beta
6//! `betai` (mirrored for negative `t`), the quantile inverts the CDF by
7//! bisection, and sampling forms `Z / sqrt(χ²_ν / ν)`.
8//!
9//! # Examples
10//!
11//! ```
12//! use stats_claw::distributions::{Cdf, Pdf};
13//! use stats_claw::distributions::TDistribution;
14//!
15//! let d = TDistribution { degrees_of_freedom: 7, ..Default::default() };
16//! // The t distribution is symmetric: density at +t equals density at -t.
17//! assert!((d.pdf(1.3) - d.pdf(-1.3)).abs() < 1e-12);
18//! // CDF at 0 is exactly 0.5.
19//! assert!((d.cdf(0.0) - 0.5).abs() < 1e-12);
20//! ```
21
22use super::super::{Cdf, LogCdf, Moments, Pdf, Quantile, Sample, bisection_quantile, count_to_f64};
23use crate::distributions::TDistribution;
24use crate::rng::SplitMix64;
25use crate::special::{betai, ln_betai_lower, ln_gamma};
26use std::f64::consts::{LN_2, PI};
27
28impl TDistribution {
29    /// Degrees of freedom as a float (`ν`).
30    fn nu(&self) -> f64 {
31        count_to_f64(self.degrees_of_freedom)
32    }
33}
34
35impl Pdf for TDistribution {
36    fn pdf(&self, x: f64) -> f64 {
37        let nu = self.nu();
38        let gamma_ratio = ln_gamma(0.5 * (nu + 1.0)) - ln_gamma(0.5 * nu);
39        let log_norm = 0.5f64.mul_add(-(nu * PI).ln(), gamma_ratio);
40        let ln_kernel = -0.5 * (nu + 1.0) * x.mul_add(x / nu, 1.0).ln();
41        (log_norm + ln_kernel).exp()
42    }
43}
44
45impl Cdf for TDistribution {
46    fn cdf(&self, x: f64) -> f64 {
47        let nu = self.nu();
48        // Iₓ with x = ν/(ν+t²) gives the two-tailed mass; halve and mirror.
49        let ib = betai(0.5 * nu, 0.5, nu / x.mul_add(x, nu));
50        if x >= 0.0 {
51            0.5f64.mul_add(-ib, 1.0)
52        } else {
53            0.5 * ib
54        }
55    }
56}
57
58impl LogCdf for TDistribution {
59    fn logsf(&self, x: f64) -> f64 {
60        // Symmetric about 0: P(T > x) = P(T < −x), so logsf(x) = logcdf(−x).
61        self.logcdf(-x)
62    }
63    fn logcdf(&self, x: f64) -> f64 {
64        let nu = self.nu();
65        // ln of the two-tailed beta mass: ln I_w(ν/2, ½) with w = ν/(ν+x²).
66        let ln_ib = ln_betai_lower(0.5 * nu, 0.5, nu / x.mul_add(x, nu));
67        if x <= 0.0 {
68            // Left tail: cdf = ½·I_w; ln cdf = ln ½ + ln I_w (finite as x → −∞).
69            ln_ib - LN_2
70        } else {
71            // Right of center: cdf = 1 − ½·I_w ≈ 1; ln_1p keeps it accurate.
72            (-0.5 * ln_ib.exp()).ln_1p()
73        }
74    }
75}
76
77impl Quantile for TDistribution {
78    fn quantile(&self, p: f64) -> f64 {
79        // The t-distribution is symmetric about 0, so anchor the median exactly
80        // and mirror the lower half onto the upper. Solving only on `[0, hi]`
81        // keeps the bracket tight and the median exact (bisection on a symmetric
82        // `[-hi, hi]` would leave an O(1e-8) residual at p = 0.5).
83        if (p - 0.5).abs() < f64::EPSILON {
84            return 0.0;
85        }
86        if p < 0.5 {
87            return -self.quantile(1.0 - p);
88        }
89        bisection_quantile(p, 0.0, 1.0e6, |x| self.cdf(x))
90    }
91}
92
93impl Moments for TDistribution {
94    fn mean(&self) -> Option<f64> {
95        if self.degrees_of_freedom > 1 {
96            Some(0.0)
97        } else {
98            None
99        }
100    }
101    fn variance(&self) -> Option<f64> {
102        let nu = self.nu();
103        if self.degrees_of_freedom > 2 {
104            Some(nu / (nu - 2.0))
105        } else {
106            None
107        }
108    }
109}
110
111impl Sample for TDistribution {
112    fn sample(&self, rng: &mut SplitMix64) -> f64 {
113        let z = rng.standard_normal();
114        let mut chi2 = 0.0;
115        for _ in 0..self.degrees_of_freedom.max(1) {
116            let g = rng.standard_normal();
117            chi2 = g.mul_add(g, chi2);
118        }
119        z / (chi2 / self.nu()).sqrt()
120    }
121}
122
123#[cfg(test)]
124mod tests {
125    use super::*;
126
127    /// The density is symmetric about zero.
128    #[test]
129    fn density_is_symmetric() {
130        let d = TDistribution {
131            degrees_of_freedom: 7,
132            ..Default::default()
133        };
134        assert!((d.pdf(1.3) - d.pdf(-1.3)).abs() < 1e-12);
135    }
136
137    /// The CDF at zero is one half.
138    #[test]
139    fn cdf_at_zero_is_half() {
140        let d = TDistribution {
141            degrees_of_freedom: 4,
142            ..Default::default()
143        };
144        assert!((d.cdf(0.0) - 0.5).abs() < 1e-12);
145    }
146
147    /// `logsf` agrees with the log of the linear survival in the body and matches
148    /// scipy `t.logsf` deep in the tail (df = 5, t = 50, where 1 − cdf underflows
149    /// past linear precision).
150    #[test]
151    fn logsf_matches_scipy_and_stays_finite() {
152        let d = TDistribution {
153            degrees_of_freedom: 5,
154            ..Default::default()
155        };
156        let body = d.logsf(1.0);
157        assert!(
158            ((body - (1.0 - d.cdf(1.0)).ln()) / body.abs()).abs() < 1e-9,
159            "logsf body was {body}"
160        );
161        let tail = d.logsf(50.0);
162        let want = -17.314_140_361_404_83; // scipy.stats.t.logsf(50, 5)
163        assert!(tail.is_finite(), "logsf(50) was {tail}");
164        assert!(
165            ((tail - want) / want).abs() < 1e-9,
166            "logsf(50) = {tail}, want {want}"
167        );
168    }
169}
stats_claw/distributions/sampling/students_t.rs

stats_claw/distributions/sampling/
students_t.rs
