qualia_core_db/solvers/learning/glm/
family.rs

1//! GLM exponential-family links — the per-family functions the IRLS loop needs.
2//! Both families use their canonical link.
3
4/// A generalized-linear-model family (canonical link).
5#[derive(Debug, Clone, Copy, PartialEq, Eq)]
6pub enum Family {
7    /// Logistic regression: Bernoulli response, logit link, `μ = σ(η)`.
8    Binomial,
9    /// Poisson regression: count response, log link, `μ = exp(η)`.
10    Poisson,
11}
12
13impl Family {
14    /// Inverse link `μ = g⁻¹(η)`.
15    pub fn inv_link(self, eta: f64) -> f64 {
16        match self {
17            // Numerically stable logistic.
18            Family::Binomial => {
19                if eta >= 0.0 {
20                    1.0 / (1.0 + (-eta).exp())
21                } else {
22                    let e = eta.exp();
23                    e / (1.0 + e)
24                }
25            }
26            Family::Poisson => eta.exp(),
27        }
28    }
29
30    /// `dμ/dη` at the current mean.
31    pub fn dmu_deta(self, mu: f64) -> f64 {
32        match self {
33            Family::Binomial => mu * (1.0 - mu),
34            Family::Poisson => mu,
35        }
36    }
37
38    /// Variance function `V(μ)`.
39    pub fn variance(self, mu: f64) -> f64 {
40        match self {
41            Family::Binomial => mu * (1.0 - mu),
42            Family::Poisson => mu,
43        }
44    }
45
46    /// Unit deviance contribution `dᵢ` (so total deviance = Σ dᵢ). Used for the
47    /// model deviance / convergence on the log-likelihood scale.
48    pub fn unit_deviance(self, y: f64, mu: f64) -> f64 {
49        const EPS: f64 = 1e-12;
50        match self {
51            Family::Binomial => {
52                let m = mu.clamp(EPS, 1.0 - EPS);
53                let a = if y > 0.0 { y * (y / m).ln() } else { 0.0 };
54                let b = if y < 1.0 {
55                    (1.0 - y) * ((1.0 - y) / (1.0 - m)).ln()
56                } else {
57                    0.0
58                };
59                2.0 * (a + b)
60            }
61            Family::Poisson => {
62                let m = mu.max(EPS);
63                let a = if y > 0.0 { y * (y / m).ln() } else { 0.0 };
64                2.0 * (a - (y - m))
65            }
66        }
67    }
68
69    /// A safe starting mean for IRLS from the raw response.
70    pub fn start_mu(self, y: f64) -> f64 {
71        match self {
72            Family::Binomial => (y + 0.5) / 2.0, // pull toward 0.5
73            Family::Poisson => (y + 0.1).max(0.1),
74        }
75    }
76}
77
78#[cfg(test)]
79mod tests {
80    use super::*;
81
82    #[test]
83    fn logistic_link_round_trip() {
84        let f = Family::Binomial;
85        assert!((f.inv_link(0.0) - 0.5).abs() < 1e-12);
86        assert!(f.inv_link(20.0) > 0.999);
87        assert!(f.inv_link(-20.0) < 0.001);
88        // variance peaks at μ=0.5.
89        assert!((f.variance(0.5) - 0.25).abs() < 1e-12);
90    }
91
92    #[test]
93    fn poisson_link() {
94        let f = Family::Poisson;
95        assert!((f.inv_link(0.0) - 1.0).abs() < 1e-12);
96        assert!((f.variance(3.0) - 3.0).abs() < 1e-12);
97    }
98
99    #[test]
100    fn deviance_is_zero_at_perfect_fit() {
101        assert!(Family::Poisson.unit_deviance(4.0, 4.0).abs() < 1e-9);
102        assert!(Family::Binomial.unit_deviance(1.0, 1.0 - 1e-13).abs() < 1e-6);
103    }
104}
qualia_core_db/solvers/learning/glm/family.rs

qualia_core_db/solvers/learning/glm/
family.rs
