stats_claw/distributions/sampling/
students_t.rs

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

stats_claw/distributions/sampling/
students_t.rs
