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}