ndarray-glm 0.1.0

//! Logistic response function where y is drawn from the Bernoulli distribution with some
//! probability p. The canonical link function is $`g(p) = \textrm{logit}(p)`$.

#[cfg(feature = "stats")]
use crate::response::Response;
use crate::{
    error::{RegressionError, RegressionResult},
    glm::{DispersionType, Glm},
    link::Link,
    math::prod_log,
    num::Float,
    response::Yval,
};
use ndarray::Array1;
#[cfg(feature = "stats")]
use statrs::distribution::Bernoulli;
use std::marker::PhantomData;

/// Logistic regression
pub struct Logistic<L = link::Logit>
where
    L: Link<Logistic<L>>,
{
    _link: PhantomData<L>,
}

/// The logistic response variable must be boolean (at least for now).
impl<L> Yval<Logistic<L>> for bool
where
    L: Link<Logistic<L>>,
{
    fn into_float<F: Float>(self) -> RegressionResult<F, F> {
        Ok(if self { F::one() } else { F::zero() })
    }
}
// Allow floats for the domain. We can't use num_traits::Float because of the
// possibility of conflicting implementations upstream, so manually implement
// for f32 and f64.
impl<L> Yval<Logistic<L>> for f32
where
    L: Link<Logistic<L>>,
{
    fn into_float<F: Float>(self) -> RegressionResult<F, F> {
        if !(0.0..=1.0).contains(&self) {
            return Err(RegressionError::InvalidY(self.to_string()));
        }
        F::from(self).ok_or_else(|| RegressionError::InvalidY(self.to_string()))
    }
}
impl<L> Yval<Logistic<L>> for f64
where
    L: Link<Logistic<L>>,
{
    fn into_float<F: Float>(self) -> RegressionResult<F, F> {
        if !(0.0..=1.0).contains(&self) {
            return Err(RegressionError::InvalidY(self.to_string()));
        }
        F::from(self).ok_or_else(|| RegressionError::InvalidY(self.to_string()))
    }
}

#[cfg(feature = "stats")]
impl<L> Response for Logistic<L>
where
    L: Link<Logistic<L>>,
{
    type DistributionType = Bernoulli;

    fn get_distribution(mu: f64, _phi: f64) -> Self::DistributionType {
        Bernoulli::new(mu).unwrap()
    }
}

/// Implementation of GLM functionality for logistic regression.
impl<L> Glm for Logistic<L>
where
    L: Link<Logistic<L>>,
{
    type Link = L;
    const DISPERSED: DispersionType = DispersionType::NoDispersion;

    /// The log of the partition function for logistic regression. The natural
    /// parameter is the logit of p.
    fn log_partition<F: Float>(nat_par: F) -> F {
        num_traits::Float::exp(nat_par).ln_1p()
    }

    /// var = mu*(1-mu)
    fn variance<F: Float>(mean: F) -> F {
        mean * (F::one() - mean)
    }

    /// This function is specialized over the default provided by Glm in order
    /// to handle over/underflow issues more precisely.
    fn log_like_natural<F>(y: F, logit_p: F) -> F
    where
        F: Float,
    {
        let (yt, xt) = if logit_p < F::zero() {
            (y, logit_p)
        } else {
            (F::one() - y, -logit_p)
        };
        yt * xt - num_traits::Float::exp(xt).ln_1p()
    }

    /// The saturated log-likelihood is zero for logistic regression when y = 0 or 1 but is less
    /// than zero for 0 < y < 1.
    fn log_like_sat<F: Float>(y: F) -> F {
        prod_log(y) + prod_log(F::one() - y)
    }
}

pub mod link {
    //! Link functions for logistic regression
    use super::*;
    use crate::link::{Canonical, Link, Transform};
    use crate::num::Float;

    /// The canonical link function for logistic regression is the logit function g(p) =
    /// log(p/(1-p)).
    pub struct Logit {}
    impl Canonical for Logit {}
    impl Link<Logistic<Logit>> for Logit {
        fn func<F: Float>(y: F) -> F {
            num_traits::Float::ln(y / (F::one() - y))
        }
        fn func_inv<F: Float>(lin_pred: F) -> F {
            (F::one() + num_traits::Float::exp(-lin_pred)).recip()
        }
    }

    /// The complementary log-log link g(p) = log(-log(1-p)) is appropriate when
    /// modeling the probability of non-zero counts when the counts are
    /// Poisson-distributed with mean lambda = exp(lin_pred).
    pub struct Cloglog {}
    impl Link<Logistic<Cloglog>> for Cloglog {
        fn func<F: Float>(y: F) -> F {
            num_traits::Float::ln(-F::ln_1p(-y))
        }
        // This quickly underflows to zero for inputs greater than ~2.
        fn func_inv<F: Float>(lin_pred: F) -> F {
            -F::exp_m1(-num_traits::Float::exp(lin_pred))
        }
    }
    impl Transform for Cloglog {
        fn nat_param<F: Float>(lin_pred: Array1<F>) -> Array1<F> {
            lin_pred.mapv(|x| num_traits::Float::ln(num_traits::Float::exp(x).exp_m1()))
        }
        fn d_nat_param<F: Float>(lin_pred: &Array1<F>) -> Array1<F> {
            let neg_exp_lin = -lin_pred.mapv(num_traits::Float::exp);
            &neg_exp_lin / &neg_exp_lin.mapv(F::exp_m1)
        }
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::{error::RegressionResult, model::ModelBuilder};
    use approx::assert_abs_diff_eq;
    use ndarray::array;

    /// A simple test where the correct value for the data is known exactly.
    #[test]
    fn log_reg() -> RegressionResult<(), f64> {
        let beta = array![0., 1.0];
        let ln2 = f64::ln(2.);
        let data_x = array![[0.], [0.], [ln2], [ln2], [ln2]];
        let data_y = array![true, false, true, true, false];
        let model = ModelBuilder::<Logistic>::data(&data_y, &data_x).build()?;
        let fit = model.fit()?;
        // dbg!(fit.n_iter);
        // NOTE: This tolerance must be higher than it would ideally be.
        // Only 2 iterations are completed, so more accuracy could presumably be achieved with a
        // lower tolerance.
        assert_abs_diff_eq!(beta, fit.result, epsilon = 0.5 * f32::EPSILON as f64);
        // let lr = fit.lr_test();
        Ok(())
    }

    #[test]
    // Confirm logit closure explicitly.
    fn logit_closure() {
        use crate::link::TestLink;
        // Because floats lose precision on difference from 1 relative to 0, higher values get
        // mapped back to infinity under closure. This is sort of fundamental to the logit
        // function and I'm not sure there's a good way around it.
        let x = array![-500., -50., -2.0, -0.2, 0., 0.5, 20.];
        super::link::Logit::check_closure(&x);
        let y = array![0., 1e-5, 0.25, 0.5, 0.8, 0.9999, 1.0];
        super::link::Logit::check_closure_y(&y);
    }

    // verify that the link and inverse are indeed inverses for the cloglog link.
    #[test]
    fn cloglog_closure() {
        use crate::link::TestLink;
        use link::Cloglog;
        let mu_test_vals = array![1e-8, 0.01, 0.1, 0.3, 0.5, 0.7, 0.9, 0.99, 0.9999999];
        Cloglog::check_closure_y(&mu_test_vals);
        let lin_test_vals = array![-10., -2., -0.1, 0.0, 0.1, 1., 2.];
        Cloglog::check_closure(&lin_test_vals);
    }

    #[test]
    fn cloglog_nat_par() {
        use crate::link::TestLink;
        use link::Cloglog;
        // nat_param(ω) = logit(cloglog^{-1}(ω)) = g_0(g^{-1}(ω))
        let lin_test_vals = array![-10., -2., -0.1, 0.0, 0.1, 1., 2.];
        Cloglog::check_nat_par::<Logistic<link::Logit>>(&lin_test_vals);
        Cloglog::check_nat_par_d(&lin_test_vals);
    }
}