ndarray_glm/response/
logistic.rs

1//! Logistic response function where y is drawn from the Bernoulli distribution with some
2//! probability p. The canonical link function is $`g(p) = \textrm{logit}(p)`$.
3
4#[cfg(feature = "stats")]
5use crate::response::Response;
6use crate::{
7    error::{RegressionError, RegressionResult},
8    glm::{DispersionType, Glm},
9    link::Link,
10    math::prod_log,
11    num::Float,
12    response::Yval,
13};
14use ndarray::Array1;
15#[cfg(feature = "stats")]
16use statrs::distribution::Bernoulli;
17use std::marker::PhantomData;
18
19/// Logistic regression
20pub struct Logistic<L = link::Logit>
21where
22    L: Link<Logistic<L>>,
23{
24    _link: PhantomData<L>,
25}
26
27/// The logistic response variable must be boolean (at least for now).
28impl<L> Yval<Logistic<L>> for bool
29where
30    L: Link<Logistic<L>>,
31{
32    fn into_float<F: Float>(self) -> RegressionResult<F, F> {
33        Ok(if self { F::one() } else { F::zero() })
34    }
35}
36// Allow floats for the domain. We can't use num_traits::Float because of the
37// possibility of conflicting implementations upstream, so manually implement
38// for f32 and f64.
39impl<L> Yval<Logistic<L>> for f32
40where
41    L: Link<Logistic<L>>,
42{
43    fn into_float<F: Float>(self) -> RegressionResult<F, F> {
44        if !(0.0..=1.0).contains(&self) {
45            return Err(RegressionError::InvalidY(self.to_string()));
46        }
47        F::from(self).ok_or_else(|| RegressionError::InvalidY(self.to_string()))
48    }
49}
50impl<L> Yval<Logistic<L>> for f64
51where
52    L: Link<Logistic<L>>,
53{
54    fn into_float<F: Float>(self) -> RegressionResult<F, F> {
55        if !(0.0..=1.0).contains(&self) {
56            return Err(RegressionError::InvalidY(self.to_string()));
57        }
58        F::from(self).ok_or_else(|| RegressionError::InvalidY(self.to_string()))
59    }
60}
61
62#[cfg(feature = "stats")]
63impl<L> Response for Logistic<L>
64where
65    L: Link<Logistic<L>>,
66{
67    type DistributionType = Bernoulli;
68
69    fn get_distribution(mu: f64, _phi: f64) -> Self::DistributionType {
70        Bernoulli::new(mu).unwrap()
71    }
72}
73
74/// Implementation of GLM functionality for logistic regression.
75impl<L> Glm for Logistic<L>
76where
77    L: Link<Logistic<L>>,
78{
79    type Link = L;
80    const DISPERSED: DispersionType = DispersionType::NoDispersion;
81
82    /// The log of the partition function for logistic regression. The natural
83    /// parameter is the logit of p.
84    fn log_partition<F: Float>(nat_par: F) -> F {
85        num_traits::Float::exp(nat_par).ln_1p()
86    }
87
88    /// var = mu*(1-mu)
89    fn variance<F: Float>(mean: F) -> F {
90        mean * (F::one() - mean)
91    }
92
93    /// This function is specialized over the default provided by Glm in order
94    /// to handle over/underflow issues more precisely.
95    fn log_like_natural<F>(y: F, logit_p: F) -> F
96    where
97        F: Float,
98    {
99        let (yt, xt) = if logit_p < F::zero() {
100            (y, logit_p)
101        } else {
102            (F::one() - y, -logit_p)
103        };
104        yt * xt - num_traits::Float::exp(xt).ln_1p()
105    }
106
107    /// The saturated log-likelihood is zero for logistic regression when y = 0 or 1 but is less
108    /// than zero for 0 < y < 1.
109    fn log_like_sat<F: Float>(y: F) -> F {
110        prod_log(y) + prod_log(F::one() - y)
111    }
112}
113
114pub mod link {
115    //! Link functions for logistic regression
116    use super::*;
117    use crate::link::{Canonical, Link, Transform};
118    use crate::num::Float;
119
120    /// The canonical link function for logistic regression is the logit function g(p) =
121    /// log(p/(1-p)).
122    pub struct Logit {}
123    impl Canonical for Logit {}
124    impl Link<Logistic<Logit>> for Logit {
125        fn func<F: Float>(y: F) -> F {
126            num_traits::Float::ln(y / (F::one() - y))
127        }
128        fn func_inv<F: Float>(lin_pred: F) -> F {
129            (F::one() + num_traits::Float::exp(-lin_pred)).recip()
130        }
131    }
132
133    /// The complementary log-log link g(p) = log(-log(1-p)) is appropriate when
134    /// modeling the probability of non-zero counts when the counts are
135    /// Poisson-distributed with mean lambda = exp(lin_pred).
136    pub struct Cloglog {}
137    impl Link<Logistic<Cloglog>> for Cloglog {
138        fn func<F: Float>(y: F) -> F {
139            num_traits::Float::ln(-F::ln_1p(-y))
140        }
141        // This quickly underflows to zero for inputs greater than ~2.
142        fn func_inv<F: Float>(lin_pred: F) -> F {
143            -F::exp_m1(-num_traits::Float::exp(lin_pred))
144        }
145    }
146    impl Transform for Cloglog {
147        fn nat_param<F: Float>(lin_pred: Array1<F>) -> Array1<F> {
148            lin_pred.mapv(|x| num_traits::Float::ln(num_traits::Float::exp(x).exp_m1()))
149        }
150        fn d_nat_param<F: Float>(lin_pred: &Array1<F>) -> Array1<F> {
151            let neg_exp_lin = -lin_pred.mapv(num_traits::Float::exp);
152            &neg_exp_lin / &neg_exp_lin.mapv(F::exp_m1)
153        }
154    }
155}
156
157#[cfg(test)]
158mod tests {
159    use super::*;
160    use crate::{error::RegressionResult, model::ModelBuilder};
161    use approx::assert_abs_diff_eq;
162    use ndarray::array;
163
164    /// A simple test where the correct value for the data is known exactly.
165    #[test]
166    fn log_reg() -> RegressionResult<(), f64> {
167        let beta = array![0., 1.0];
168        let ln2 = f64::ln(2.);
169        let data_x = array![[0.], [0.], [ln2], [ln2], [ln2]];
170        let data_y = array![true, false, true, true, false];
171        let model = ModelBuilder::<Logistic>::data(&data_y, &data_x).build()?;
172        let fit = model.fit()?;
173        // dbg!(fit.n_iter);
174        // NOTE: This tolerance must be higher than it would ideally be.
175        // Only 2 iterations are completed, so more accuracy could presumably be achieved with a
176        // lower tolerance.
177        assert_abs_diff_eq!(beta, fit.result, epsilon = 0.5 * f32::EPSILON as f64);
178        // let lr = fit.lr_test();
179        Ok(())
180    }
181
182    #[test]
183    // Confirm logit closure explicitly.
184    fn logit_closure() {
185        use crate::link::TestLink;
186        // Because floats lose precision on difference from 1 relative to 0, higher values get
187        // mapped back to infinity under closure. This is sort of fundamental to the logit
188        // function and I'm not sure there's a good way around it.
189        let x = array![-500., -50., -2.0, -0.2, 0., 0.5, 20.];
190        super::link::Logit::check_closure(&x);
191        let y = array![0., 1e-5, 0.25, 0.5, 0.8, 0.9999, 1.0];
192        super::link::Logit::check_closure_y(&y);
193    }
194
195    // verify that the link and inverse are indeed inverses for the cloglog link.
196    #[test]
197    fn cloglog_closure() {
198        use crate::link::TestLink;
199        use link::Cloglog;
200        let mu_test_vals = array![1e-8, 0.01, 0.1, 0.3, 0.5, 0.7, 0.9, 0.99, 0.9999999];
201        Cloglog::check_closure_y(&mu_test_vals);
202        let lin_test_vals = array![-10., -2., -0.1, 0.0, 0.1, 1., 2.];
203        Cloglog::check_closure(&lin_test_vals);
204    }
205
206    #[test]
207    fn cloglog_nat_par() {
208        use crate::link::TestLink;
209        use link::Cloglog;
210        // nat_param(ω) = logit(cloglog^{-1}(ω)) = g_0(g^{-1}(ω))
211        let lin_test_vals = array![-10., -2., -0.1, 0.0, 0.1, 1., 2.];
212        Cloglog::check_nat_par::<Logistic<link::Logit>>(&lin_test_vals);
213        Cloglog::check_nat_par_d(&lin_test_vals);
214    }
215}
ndarray_glm/response/logistic.rs

ndarray_glm/response/
logistic.rs
