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

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

/// Implements the [Student's T](https://en.wikipedia.org/wiki/Student%27s_t-distribution) distribution.
#[derive(Debug, Clone, Copy)]
pub struct T {
    /// Degrees of freedom
    dof: f64,
}

pub type StudentsT = T;

impl T {
    /// Create a new t distribution with
    ///
    /// # Errors
    /// Panics if degrees of freedom is not positive.
    pub fn new(dof: f64) -> Self {
        assert!(dof > 0., "Degrees of freedom must be positive.");
        T { dof }
    }
    pub fn set_dof(&mut self, dof: f64) -> &mut Self {
        assert!(dof > 0., "Degrees of freedom must be positive.");
        self.dof = dof;
        self
    }
}

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

impl Distribution for T {
    type Output = f64;
    /// Samples from the given T distribution.
    fn sample(&self) -> f64 {
        (self.dof / 2.).sqrt() * Normal::default().sample()
            / Gamma::new(self.dof / 2., 1.).sample().sqrt()
    }
}

impl Distribution1D for T {
    fn update(&mut self, params: &[f64]) {
        self.set_dof(params[0]);
    }
}

impl Continuous for T {
	type PDFType = f64;
    /// Calculates the probability density function for the given T distribution at `x`.
    fn pdf(&self, x: f64) -> f64 {
        gamma((self.dof + 1.) / 2.)
            / ((self.dof * std::f64::consts::PI).sqrt() * gamma(self.dof / 2.))
            * (1. + x.powi(2) / self.dof).powf(-(self.dof - 1.) / 2.)
    }
}

impl Mean for T {
    type MeanType = f64;
    /// Calculates the mean of the T distribution, which is 0 when the degrees of freedom is
    /// greater than 1, and undefined otherwise.
    ///
    fn mean(&self) -> f64 {
        if self.dof > 1. {
            0.
        } else {
            f64::NAN
        }
    }
}

impl Variance for T {
    type VarianceType = f64;
    /// Calculates the variance of the T distribution.
    ///
    /// # Remarks
    /// This is not defined when degrees of freedom is less than or equal to 1, and infinity when
    /// degrees of freedom is in (1, 2].
    fn var(&self) -> f64 {
        if self.dof > 2. {
            self.dof / (self.dof - 2.)
        } else if (1. < self.dof) & (self.dof <= 2.) {
            f64::INFINITY
        } else {
            f64::NAN
        }
    }
}

#[cfg(test)]
mod tests {

    use super::*;
    use crate::statistics::mean;
    use approx_eq::assert_approx_eq;

    #[test]
    fn test_moments() {
        let t = T::new(2.);
        let data = t.sample_n(1e6 as usize);
        assert_approx_eq!(mean(&data), 0., 1e-2);
    }
}