r2rs-stats 0.1.1

Statistics programming for Rust based on R's stats package
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81

// "Whatever you do, work at it with all your heart, as working for the Lord,
// not for human masters, since you know that you will receive an inheritance
// from the Lord as a reward. It is the Lord Christ you are serving."
// (Col 3:23-24)

use nalgebra::DMatrix;
use num_traits::Float;

use crate::regression::glm::link::{core::Link, identity::Identity, log::Log};

#[derive(Clone, Debug)]
pub struct Power {
    lambda: f64,
}

impl Default for Power {
    fn default() -> Power {
        Self { lambda: 1.0 }.with_lambda(1.0)
    }
}

impl Power {
    pub fn new() -> Self {
        Self::default()
    }

    pub fn with_lambda(self, lambda: f64) -> Self {
        Self { lambda }
    }
}

impl Link for Power {
    fn link_func(&self, mu: &DMatrix<f64>) -> DMatrix<f64> {
        if self.lambda <= 0.0 {
            Log::new().link_func(mu)
        } else if self.lambda == 1.0 {
            Identity::new().link_func(mu)
        } else {
            let mut ret = mu.clone();
            ret.iter_mut().for_each(|mu| *mu = mu.powf(self.lambda));
            ret
        }
    }

    fn link_inverse(&self, eta: &DMatrix<f64>) -> DMatrix<f64> {
        if self.lambda <= 0.0 {
            Log::new().link_inverse(eta)
        } else if self.lambda == 1.0 {
            Identity::new().link_inverse(eta)
        } else {
            let mut ret = eta.clone();
            ret.iter_mut()
                .for_each(|eta| *eta = eta.powf(1.0 / self.lambda).max(f64::epsilon()));
            ret
        }
    }

    fn mu_eta(&self, eta: &DMatrix<f64>) -> DMatrix<f64> {
        if self.lambda <= 0.0 {
            Log::new().mu_eta(eta)
        } else if self.lambda == 1.0 {
            Identity::new().mu_eta(eta)
        } else {
            let mut ret = eta.clone();
            ret.iter_mut().for_each(|eta| {
                *eta = ((1.0 / self.lambda) * eta.powf(1.0 / self.lambda - 1.0)).max(f64::epsilon())
            });
            ret
        }
    }

    fn valid_eta(&self, eta: &DMatrix<f64>) -> bool {
        if self.lambda <= 0.0 {
            Log::new().valid_eta(eta)
        } else if self.lambda == 1.0 {
            Identity::new().valid_eta(eta)
        } else {
            eta.iter().all(|&eta| eta.is_finite() && eta > 0.0)
        }
    }
}