qualia_core_db/solvers/statistics/distributions/
fisher_f.rs

1//! Fisher–Snedecor F-distribution — pdf / cdf / quantile + upper-tail p-value, used
2//! by ANOVA and variance-ratio tests. CDF is exact via the regularized incomplete
3//! beta `I_x(d₁/2, d₂/2)` with `x = d₁f/(d₁f + d₂)`.
4
5use super::special::{betai, ln_gamma};
6
7/// pdf with numerator/denominator dof `d1, d2 > 0`, `x ≥ 0`.
8pub fn pdf(x: f64, d1: f64, d2: f64) -> f64 {
9    debug_assert!(d1 > 0.0 && d2 > 0.0);
10    if x <= 0.0 {
11        return 0.0;
12    }
13    // ln of (d1/d2)^{d1/2} x^{d1/2-1} (1+d1 x/d2)^{-(d1+d2)/2} / B(d1/2,d2/2)
14    let ln_b = ln_gamma(d1 / 2.0) + ln_gamma(d2 / 2.0) - ln_gamma((d1 + d2) / 2.0);
15    let ln_num = (d1 / 2.0) * (d1 / d2).ln() + (d1 / 2.0 - 1.0) * x.ln()
16        - ((d1 + d2) / 2.0) * (1.0 + d1 * x / d2).ln();
17    (ln_num - ln_b).exp()
18}
19
20/// cdf `P(X ≤ x)` = `I_x(d₁/2, d₂/2)`, `x = d₁f/(d₁f + d₂)`.
21pub fn cdf(f: f64, d1: f64, d2: f64) -> f64 {
22    debug_assert!(d1 > 0.0 && d2 > 0.0);
23    if f <= 0.0 {
24        return 0.0;
25    }
26    let x = d1 * f / (d1 * f + d2);
27    betai(d1 / 2.0, d2 / 2.0, x)
28}
29
30/// Upper-tail p-value `P(X ≥ f)` — the ANOVA / variance-ratio p-value.
31pub fn upper_p(f: f64, d1: f64, d2: f64) -> f64 {
32    1.0 - cdf(f, d1, d2)
33}
34
35/// Inverse cdf (quantile) for `0 < p < 1`.
36pub fn quantile(p: f64, d1: f64, d2: f64) -> f64 {
37    super::invert_cdf(p, Some(0.0), |f| cdf(f, d1, d2))
38}
39
40#[cfg(test)]
41mod tests {
42    use super::*;
43
44    #[test]
45    fn known_critical_values() {
46        // F_{0.95}(5, 10) = 3.3258; F_{0.95}(1, 10) = 4.9646.
47        assert!((quantile(0.95, 5.0, 10.0) - 3.325_835).abs() < 1e-3);
48        assert!((quantile(0.95, 1.0, 10.0) - 4.964_603).abs() < 1e-3);
49        // Upper-tail p of that critical value ≈ 0.05.
50        assert!((upper_p(3.325_835, 5.0, 10.0) - 0.05).abs() < 1e-4);
51    }
52
53    #[test]
54    fn cdf_monotone_and_bounds() {
55        assert_eq!(cdf(0.0, 3.0, 8.0), 0.0);
56        assert!(cdf(1.0, 3.0, 8.0) < cdf(5.0, 3.0, 8.0));
57        assert!(cdf(1e6, 3.0, 8.0) > 0.999);
58    }
59
60    #[test]
61    fn quantile_inverts_cdf() {
62        for &p in &[0.1, 0.5, 0.9, 0.99] {
63            let f = quantile(p, 4.0, 20.0);
64            assert!((cdf(f, 4.0, 20.0) - p).abs() < 1e-7, "p={p}");
65        }
66    }
67}
qualia_core_db/solvers/statistics/distributions/fisher_f.rs

qualia_core_db/solvers/statistics/distributions/
fisher_f.rs
