compute 0.2.3

A crate for statistical computing.
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
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131

use crate::distributions::*;
use crate::functions::gamma;

/// Implements the [Gamma](https://en.wikipedia.org/wiki/Gamma_distribution) distribution.
#[derive(Debug, Clone, Copy)]
pub struct Gamma {
    /// Shape parameter α.
    alpha: f64,
    /// Rate parameter β.
    beta: f64,
    normal_gen: Normal,
    uniform_gen: Uniform,
}

impl Gamma {
    /// Create a new Gamma distribution with shape `alpha` and rate `beta`.
    ///
    /// # Errors
    /// Panics if `alpha <= 0` or `beta <= 0`.
    pub fn new(alpha: f64, beta: f64) -> Self {
        if alpha <= 0. || beta <= 0. {
            panic!("Both alpha and beta must be positive.");
        }
        Gamma {
            alpha,
            beta,
            normal_gen: Normal::new(0., 1.),
            uniform_gen: Uniform::new(0., 1.),
        }
    }
    pub fn set_alpha(&mut self, alpha: f64) -> &mut Self {
        if alpha <= 0. {
            panic!("Alpha must be positive.");
        }
        self.alpha = alpha;
        self
    }
    pub fn set_beta(&mut self, beta: f64) -> &mut Self {
        if beta <= 0. {
            panic!("Beta must be positive.");
        }
        self.beta = beta;
        self
    }
}

impl Default for Gamma {
    fn default() -> Self {
        Self::new(1., 1.)
    }
}

impl Distribution for Gamma {
    type Output = f64;
    /// Samples from the given Gamma distribution.
    ///
    /// # Remarks
    /// Uses the algorithm from Marsaglia and Tsang 2000. Applies the squeeze
    /// method and has nearly constant average time for `alpha >= 1`.
    fn sample(&self) -> f64 {
        let d = self.alpha - 1. / 3.;
        loop {
            let (x, v) = loop {
                let x = self.normal_gen.sample();
                let v = (1. + x / (9. * d).sqrt()).powi(3);
                if v > 0. {
                    break (x, v);
                }
            };
            let u = self.uniform_gen.sample();
            if u < 1. - 0.0331 * x.powi(4) {
                return d * v / self.beta;
            }
            if u.ln() < 0.5 * x.powi(2) + d * (1. - v + v.ln()) {
                return d * v / self.beta;
            }
        }
    }
}

impl Distribution1D for Gamma {
    fn update(&mut self, params: &[f64]) {
        self.set_alpha(params[0]).set_beta(params[1]);
    }
}

impl Continuous for Gamma {
	type PDFType = f64;
    /// Calculates the probability density function for the given Gamma function at `x`.
    ///
    /// # Remarks
    /// x should be positive.
    fn pdf(&self, x: f64) -> f64 {
        if x <= 0. {
            return 0.;
        }
        self.beta.powf(self.alpha) / gamma(self.alpha)
            * x.powf(self.alpha - 1.)
            * (-self.beta * x).exp()
    }
}

impl Mean for Gamma {
    type MeanType = f64;
    /// Calculates the mean, which for a Gamma(a, b) distribution is given by `a / b`.
    fn mean(&self) -> f64 {
        self.alpha / self.beta
    }
}

impl Variance for Gamma {
    type VarianceType = f64;
    /// Calculates the variance of the given Gamma distribution.
    fn var(&self) -> f64 {
        self.alpha / self.beta.powi(2)
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::statistics::{mean, var};
    use approx_eq::assert_approx_eq;

    #[test]
    fn test_moments() {
        let data = Gamma::new(2., 4.).sample_n(1e6 as usize);
        assert_approx_eq!(0.5, mean(&data), 1e-2);
        assert_approx_eq!(0.125, var(&data), 1e-2);
    }
}