numra_stats/distributions/
f_dist.rs

1//! F-distribution.
2//!
3//! Author: Moussa Leblouba
4//! Date: 9 February 2026
5//! Modified: 2 May 2026
6
7use numra_core::Scalar;
8use numra_special::{betainc, lgamma};
9use rand::RngCore;
10
11use super::gamma_dist::GammaDist;
12use super::ContinuousDistribution;
13
14/// F-distribution with d1 and d2 degrees of freedom.
15#[derive(Clone, Debug)]
16pub struct FDist<S: Scalar> {
17    pub df1: S,
18    pub df2: S,
19}
20
21impl<S: Scalar> FDist<S> {
22    pub fn new(df1: S, df2: S) -> Self {
23        Self { df1, df2 }
24    }
25}
26
27impl<S: Scalar> ContinuousDistribution<S> for FDist<S> {
28    fn pdf(&self, x: S) -> S {
29        if x <= S::ZERO {
30            return S::ZERO;
31        }
32        let half = S::HALF;
33        let d1 = self.df1;
34        let d2 = self.df2;
35        let log_pdf = half * d1 * (d1 * x / (d1 * x + d2)).ln()
36            + half * d2 * (d2 / (d1 * x + d2)).ln()
37            - x.ln()
38            - lbeta(d1 * half, d2 * half);
39        log_pdf.exp()
40    }
41
42    fn cdf(&self, x: S) -> S {
43        if x <= S::ZERO {
44            return S::ZERO;
45        }
46        let half = S::HALF;
47        let d1 = self.df1;
48        let d2 = self.df2;
49        let z = d1 * x / (d1 * x + d2);
50        betainc(d1 * half, d2 * half, z)
51    }
52
53    fn quantile(&self, p: S) -> S {
54        if p <= S::ZERO {
55            return S::ZERO;
56        }
57        if p >= S::ONE {
58            return S::INFINITY;
59        }
60        // Newton iteration
61        let mut x = S::ONE; // initial guess
62        for _ in 0..100 {
63            let f_val = self.cdf(x) - p;
64            let f_prime = self.pdf(x);
65            if f_prime.to_f64().abs() < 1e-300 {
66                break;
67            }
68            let step = f_val / f_prime;
69            x -= step;
70            if x <= S::ZERO {
71                x = S::from_f64(1e-10);
72            }
73            if step.to_f64().abs() < 1e-12 * x.to_f64().abs() {
74                break;
75            }
76        }
77        x
78    }
79
80    fn mean(&self) -> S {
81        let two = S::TWO;
82        if self.df2 > two {
83            self.df2 / (self.df2 - two)
84        } else {
85            S::INFINITY
86        }
87    }
88
89    fn variance(&self) -> S {
90        let two = S::TWO;
91        let four = S::from_f64(4.0);
92        if self.df2 > four {
93            let d1 = self.df1;
94            let d2 = self.df2;
95            two * d2 * d2 * (d1 + d2 - two) / (d1 * (d2 - two) * (d2 - two) * (d2 - four))
96        } else {
97            S::INFINITY
98        }
99    }
100
101    fn sample(&self, rng: &mut dyn RngCore) -> S {
102        // F = (X1/d1) / (X2/d2) where Xi ~ Chi2(di)
103        let half = S::HALF;
104        let g1 = GammaDist::new(self.df1 * half, half);
105        let g2 = GammaDist::new(self.df2 * half, half);
106        let x1 = g1.sample(rng);
107        let x2 = g2.sample(rng);
108        (x1 / self.df1) / (x2 / self.df2)
109    }
110}
111
112fn lbeta<S: Scalar>(a: S, b: S) -> S {
113    lgamma(a) + lgamma(b) - lgamma(a + b)
114}
115
116#[cfg(test)]
117mod tests {
118    use super::*;
119
120    #[test]
121    fn test_f_dist_mean() {
122        let f = FDist::new(5.0_f64, 10.0);
123        assert!((f.mean() - 10.0 / 8.0).abs() < 1e-12);
124    }
125
126    #[test]
127    fn test_f_dist_cdf_at_zero() {
128        let f = FDist::new(3.0_f64, 5.0);
129        assert!(f.cdf(0.0).abs() < 1e-10);
130    }
131
132    #[test]
133    fn test_f_dist_quantile_roundtrip() {
134        let f = FDist::new(5.0_f64, 10.0);
135        for &p in &[0.1, 0.5, 0.9] {
136            let x = f.quantile(p);
137            let p2 = f.cdf(x);
138            assert!((p - p2).abs() < 1e-5, "p={}, p2={}", p, p2);
139        }
140    }
141}
numra_stats/distributions/f_dist.rs

numra_stats/distributions/
f_dist.rs
