prospicio_models/
link.rs

1//! Link functions: `η = g(μ)` between the mean and the linear predictor.
2
3use prospicio_math::special::{norm_cdf, norm_pdf, norm_quantile};
4
5/// A link function `η = g(μ)`, with its inverse and derivative.
6///
7/// A closed enum, so a fitted model's link serializes as data.
8///
9/// ```
10/// use prospicio_models::Link;
11///
12/// let g = Link::Log;
13/// assert!((g.inverse(g.link(3.0)) - 3.0).abs() < 1e-15);
14/// assert_eq!(g.mu_eta(0.0), 1.0);
15/// ```
16#[derive(Debug, Clone, Copy, PartialEq)]
17pub enum Link {
18    /// `η = μ`.
19    Identity,
20    /// `η = ln μ`: multiplicative models, the actuarial default.
21    Log,
22    /// `η = ln(μ / (1 - μ))`.
23    Logit,
24    /// `η = Φ⁻¹(μ)`.
25    Probit,
26    /// `η = ln(-ln(1 - μ))`.
27    Cloglog,
28    /// `η = 1 / μ`, the gamma family's canonical link.
29    Inverse,
30    /// `η = 1 / μ²`, the inverse Gaussian's canonical link.
31    InverseSquared,
32    /// `η = μ^λ` for `λ ≠ 0` (`λ = 0` is the log link).
33    Power(f64),
34}
35
36impl Link {
37    /// `η = g(μ)`.
38    pub fn link(&self, mu: f64) -> f64 {
39        match *self {
40            Self::Identity => mu,
41            Self::Log => mu.ln(),
42            Self::Logit => (mu / (1.0 - mu)).ln(),
43            Self::Probit => norm_quantile(mu),
44            Self::Cloglog => (-(-mu).ln_1p()).ln(),
45            Self::Inverse => 1.0 / mu,
46            Self::InverseSquared => 1.0 / (mu * mu),
47            Self::Power(0.0) => mu.ln(),
48            Self::Power(l) => mu.powf(l),
49        }
50    }
51
52    /// `μ = g⁻¹(η)`.
53    pub fn inverse(&self, eta: f64) -> f64 {
54        match *self {
55            Self::Identity => eta,
56            Self::Log => eta.exp(),
57            Self::Logit => 1.0 / (1.0 + (-eta).exp()),
58            Self::Probit => norm_cdf(eta),
59            Self::Cloglog => -(-eta.exp()).exp_m1(),
60            Self::Inverse => 1.0 / eta,
61            Self::InverseSquared => 1.0 / eta.sqrt(),
62            Self::Power(0.0) => eta.exp(),
63            Self::Power(l) => eta.powf(1.0 / l),
64        }
65    }
66
67    /// `dμ/dη` at `η`.
68    pub fn mu_eta(&self, eta: f64) -> f64 {
69        match *self {
70            Self::Identity => 1.0,
71            Self::Log => eta.exp(),
72            Self::Logit => {
73                let e = (-eta.abs()).exp();
74                e / ((1.0 + e) * (1.0 + e))
75            }
76            Self::Probit => norm_pdf(eta),
77            Self::Cloglog => (eta - eta.exp()).exp(),
78            Self::Inverse => -1.0 / (eta * eta),
79            Self::InverseSquared => -0.5 * eta.powf(-1.5),
80            Self::Power(0.0) => eta.exp(),
81            Self::Power(l) => eta.powf(1.0 / l - 1.0) / l,
82        }
83    }
84
85    /// `d²μ/dη²` at `η`, for the observed information of a non-canonical
86    /// link.
87    pub fn mu_eta2(&self, eta: f64) -> f64 {
88        match *self {
89            Self::Identity => 0.0,
90            Self::Log | Self::Power(0.0) => eta.exp(),
91            Self::Logit => {
92                let mu = self.inverse(eta);
93                self.mu_eta(eta) * (1.0 - 2.0 * mu)
94            }
95            Self::Probit => -eta * norm_pdf(eta),
96            Self::Cloglog => self.mu_eta(eta) * (1.0 - eta.exp()),
97            Self::Inverse => 2.0 / (eta * eta * eta),
98            Self::InverseSquared => 0.75 * eta.powf(-2.5),
99            Self::Power(l) => (1.0 / l - 1.0) / l * eta.powf(1.0 / l - 2.0),
100        }
101    }
102
103    /// Whether `η` maps to a mean the link can produce (for example, a
104    /// positive `η` under the inverse link).
105    pub fn valid_eta(&self, eta: f64) -> bool {
106        eta.is_finite()
107            && match *self {
108                Self::Inverse | Self::InverseSquared => eta > 0.0,
109                Self::Power(l) if l != 0.0 => eta > 0.0,
110                _ => true,
111            }
112    }
113}
114
115#[cfg(test)]
116mod tests {
117    use super::*;
118
119    #[test]
120    fn inverses_and_derivatives_agree() {
121        let links = [
122            (Link::Identity, 2.5),
123            (Link::Log, 2.5),
124            (Link::Logit, 0.3),
125            (Link::Probit, 0.3),
126            (Link::Cloglog, 0.3),
127            (Link::Inverse, 2.5),
128            (Link::InverseSquared, 2.5),
129            (Link::Power(0.5), 2.5),
130            (Link::Power(0.0), 2.5),
131        ];
132        for (g, mu) in links {
133            let eta = g.link(mu);
134            assert!((g.inverse(eta) - mu).abs() < 1e-13, "{g:?}");
135            let h = 1e-6;
136            let numeric = (g.inverse(eta + h) - g.inverse(eta - h)) / (2.0 * h);
137            assert!(
138                (g.mu_eta(eta) - numeric).abs() < 1e-7 * numeric.abs().max(1.0),
139                "{g:?}"
140            );
141            let numeric2 = (g.mu_eta(eta + h) - g.mu_eta(eta - h)) / (2.0 * h);
142            assert!(
143                (g.mu_eta2(eta) - numeric2).abs() < 1e-7 * numeric2.abs().max(1.0),
144                "{g:?}"
145            );
146        }
147    }
148}
prospicio_models/link.rs

prospicio_models/
link.rs
