ndarray_glm/
link.rs

1//! Defines traits for link functions
2
3use crate::{glm::Glm, num::Float};
4use ndarray::Array1;
5
6/// Describes the functions to map to and from the linear predictors and the
7/// expectation of the response. It is constrained mathematically by the
8/// response distribution and the transformation of the linear predictor.
9// TODO: The link function and its inverse are independent of the response
10// distribution. This could be refactored to separate the function itself from
11// the transformation that works with the distribution.
12pub trait Link<M: Glm>: Transform {
13    /// Maps the expectation value of the response variable to the linear
14    /// predictor. In general this is determined by a composition of the inverse
15    /// natural parameter transformation and the canonical link function.
16    fn func<F: Float>(y: F) -> F;
17    // fn func<F: Float>(y: Array1<F>) -> Array1<F>;
18    /// Maps the linear predictor to the expectation value of the response.
19    // TODO: There may not be a point in using Array versions of these functions
20    // since clones are necessary anyway. Perhaps we could simply define the
21    // scalar function and use mapv().
22    fn func_inv<F: Float>(lin_pred: F) -> F;
23    // fn func_inv<F: Float>(lin_pred: Array1<F>) -> Array1<F>;
24}
25
26pub trait Transform {
27    /// The natural parameter(s) of the response distribution as a function
28    /// of the linear predictor. For canonical link functions this is the
29    /// identity. It must be monotonic, invertible, and twice-differentiable.
30    /// For link function g and canonical link function g_0 it is equal to
31    /// g_0 ( g^{-1}(lin_pred) ) .
32    fn nat_param<F: Float>(lin_pred: Array1<F>) -> Array1<F>;
33    /// The derivative of the transformation to the natural parameter. If it is
34    /// zero in a region that the IRLS is in the algorithm may have difficulty
35    /// converging.
36    fn d_nat_param<F: Float>(lin_pred: &Array1<F>) -> Array1<F>;
37    /// Adjust the error and variance terms of the likelihood function based on
38    /// the first and second derivatives of the transformation. The adjustment
39    /// is performed simultaneously. The linear predictor must be
40    /// un-transformed, i.e. it must be X*beta without the transformation
41    /// applied.
42    fn adjust_errors_variance<F: Float>(
43        errors: Array1<F>,
44        variance: Array1<F>,
45        lin_pred: &Array1<F>,
46    ) -> (Array1<F>, Array1<F>) {
47        let eta_d = Self::d_nat_param(lin_pred);
48        let err_adj = &eta_d * &errors;
49        // The second-derivative term in the variance matrix can lead it to not
50        // be positive-definite. In fact, the second term should vanish when
51        // taking the expecation of Y to give the Fisher information.
52        // let var_adj = &eta_d * &variance * eta_d - eta_dd * errors;
53        let var_adj = &eta_d * &variance * eta_d;
54        (err_adj, var_adj)
55    }
56}
57
58/// The canonical transformation by definition equates the linear predictor with
59/// the natural parameter of the response distribution. Implementing this trait
60/// for a link function automatically defines the trivial transformation
61/// functions.
62pub trait Canonical {}
63impl<T> Transform for T
64where
65    T: Canonical,
66{
67    /// By defintion this function is the identity function for canonical links.
68    #[inline]
69    fn nat_param<F: Float>(lin_pred: Array1<F>) -> Array1<F> {
70        lin_pred
71    }
72    #[inline]
73    fn d_nat_param<F: Float>(lin_pred: &Array1<F>) -> Array1<F> {
74        Array1::<F>::ones(lin_pred.len())
75    }
76    /// The canonical link function requires no transformation of the error and variance terms.
77    #[inline]
78    fn adjust_errors_variance<F: Float>(
79        errors: Array1<F>,
80        variance: Array1<F>,
81        _lin_pred: &Array1<F>,
82    ) -> (Array1<F>, Array1<F>) {
83        (errors, variance)
84    }
85}
ndarray_glm/link.rs

ndarray_glm/
link.rs
