stats_claw/distributions/sampling/
f_dist.rs

1//! F-distribution numerics, for the [`FDistribution`] parameter struct.
2//!
3//! Equivalent to `scipy.stats.f(numerator_df, denominator_df)`. The density uses
4//! `ln_beta`, the CDF is the regularized incomplete beta `betai`, the quantile
5//! inverts it by bisection, and sampling forms the ratio of two scaled
6//! chi-squared draws, `(χ²_{d1}/d1) / (χ²_{d2}/d2)`.
7//!
8//! # Examples
9//!
10//! ```
11//! use stats_claw::distributions::{Cdf, Moments};
12//! use stats_claw::distributions::FDistribution;
13//!
14//! let d = FDistribution { numerator_df: 6, denominator_df: 12, ..Default::default() };
15//! // Mean is d2 / (d2 - 2) = 12 / 10 = 1.2.
16//! assert!(matches!(d.mean(), Some(m) if (m - 1.2).abs() < 1e-12));
17//! // CDF saturates toward 1 in the right tail.
18//! assert!(d.cdf(1_000.0) > 1.0 - 1e-9);
19//! ```
20
21use super::super::{Cdf, LogCdf, Moments, Pdf, Quantile, Sample, bisection_quantile, count_to_f64};
22use crate::distributions::FDistribution;
23use crate::rng::SplitMix64;
24use crate::special::{betai, ln_beta, ln_betai_lower, ln_betai_upper};
25
26impl FDistribution {
27    /// Numerator degrees of freedom as a float (`d1`).
28    fn d1(&self) -> f64 {
29        count_to_f64(self.numerator_df)
30    }
31    /// Denominator degrees of freedom as a float (`d2`).
32    fn d2(&self) -> f64 {
33        count_to_f64(self.denominator_df)
34    }
35}
36
37impl Pdf for FDistribution {
38    fn pdf(&self, x: f64) -> f64 {
39        if x <= 0.0 {
40            return 0.0;
41        }
42        let (d1, d2) = (self.d1(), self.d2());
43        // log f(x) = ½d1·ln(d1·x) + ½d2·ln d2 − ½(d1+d2)·ln(d1·x+d2) − ln B − ln x
44        let log_num = (0.5 * d1).mul_add((d1 * x).ln(), 0.5 * d2 * d2.ln());
45        let log_den =
46            (0.5 * (d1 + d2)).mul_add(d1.mul_add(x, d2).ln(), ln_beta(0.5 * d1, 0.5 * d2));
47        (log_num - log_den - x.ln()).exp()
48    }
49}
50
51impl Cdf for FDistribution {
52    fn cdf(&self, x: f64) -> f64 {
53        if x <= 0.0 {
54            return 0.0;
55        }
56        let (d1, d2) = (self.d1(), self.d2());
57        let d1x = d1 * x;
58        betai(0.5 * d1, 0.5 * d2, d1x / (d1x + d2))
59    }
60}
61
62impl LogCdf for FDistribution {
63    fn logsf(&self, x: f64) -> f64 {
64        if x <= 0.0 {
65            return 0.0;
66        }
67        let (d1, d2) = (self.d1(), self.d2());
68        let d1x = d1 * x;
69        // sf = 1 − I_w(d1/2, d2/2); the upper log tail keeps the small side fast.
70        ln_betai_upper(0.5 * d1, 0.5 * d2, d1x / (d1x + d2))
71    }
72    fn logcdf(&self, x: f64) -> f64 {
73        if x <= 0.0 {
74            return f64::NEG_INFINITY;
75        }
76        let (d1, d2) = (self.d1(), self.d2());
77        let d1x = d1 * x;
78        ln_betai_lower(0.5 * d1, 0.5 * d2, d1x / (d1x + d2))
79    }
80}
81
82impl Quantile for FDistribution {
83    fn quantile(&self, p: f64) -> f64 {
84        // Heavy right tail; a wide bracket covers the near-1 grid point.
85        bisection_quantile(p, 0.0, 1.0e7, |x| self.cdf(x))
86    }
87}
88
89impl Moments for FDistribution {
90    fn mean(&self) -> Option<f64> {
91        let d2 = self.d2();
92        if self.denominator_df > 2 {
93            Some(d2 / (d2 - 2.0))
94        } else {
95            None
96        }
97    }
98    fn variance(&self) -> Option<f64> {
99        if self.denominator_df <= 4 {
100            return None;
101        }
102        let (d1, d2) = (self.d1(), self.d2());
103        let numerator = 2.0 * d2 * d2 * (d1 + d2 - 2.0);
104        let denominator = d1 * (d2 - 2.0) * (d2 - 2.0) * (d2 - 4.0);
105        Some(numerator / denominator)
106    }
107}
108
109impl Sample for FDistribution {
110    fn sample(&self, rng: &mut SplitMix64) -> f64 {
111        let num = chi_squared_over_df(self.numerator_df, rng);
112        let den = chi_squared_over_df(self.denominator_df, rng);
113        num / den
114    }
115}
116
117/// Draws `χ²_df / df` — a scaled chi-squared with `df` degrees of freedom.
118fn chi_squared_over_df(df: i64, rng: &mut SplitMix64) -> f64 {
119    let mut sum = 0.0;
120    for _ in 0..df.max(1) {
121        let z = rng.standard_normal();
122        sum = z.mul_add(z, sum);
123    }
124    sum / count_to_f64(df.max(1))
125}
126
127#[cfg(test)]
128mod tests {
129    use super::*;
130
131    /// The mean is `d2 / (d2 − 2)` when the denominator df exceeds 2.
132    #[test]
133    fn mean_for_large_denominator() {
134        let d = FDistribution {
135            numerator_df: 6,
136            denominator_df: 12,
137            ..Default::default()
138        };
139        assert!(matches!(d.mean(), Some(m) if (m - 12.0 / 10.0).abs() < 1e-12));
140    }
141
142    /// The CDF saturates toward 1 in the upper tail.
143    #[test]
144    fn cdf_saturates() {
145        let d = FDistribution {
146            numerator_df: 6,
147            denominator_df: 12,
148            ..Default::default()
149        };
150        assert!(d.cdf(1e6) > 1.0 - 1e-9);
151    }
152
153    /// `logsf` agrees with the log of the linear survival in the body and matches
154    /// scipy `f.logsf` deep in the tail (d1 = 3, d2 = 10, x = 1e4).
155    #[test]
156    fn logsf_matches_scipy_and_stays_finite() {
157        let d = FDistribution {
158            numerator_df: 3,
159            denominator_df: 10,
160            ..Default::default()
161        };
162        let body = d.logsf(1.0);
163        assert!(
164            ((body - (1.0 - d.cdf(1.0)).ln()) / body.abs()).abs() < 1e-9,
165            "logsf body was {body}"
166        );
167        let tail = d.logsf(1e4);
168        let want = -39.037_790_534_655_89; // scipy.stats.f.logsf(1e4, 3, 10)
169        assert!(tail.is_finite(), "logsf(1e4) was {tail}");
170        assert!(
171            ((tail - want) / want).abs() < 1e-9,
172            "logsf(1e4) = {tail}, want {want}"
173        );
174    }
175}
stats_claw/distributions/sampling/f_dist.rs

stats_claw/distributions/sampling/
f_dist.rs
