statslicer 0.10.0

Guacamole provides a linearly-seekable random number generator.
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
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155

//! Statistical moments.
// TODO(rescrv):  This was copied from biometrics.  Decide how to proceed.

/// Moments are the statistical moments of mean (m1), standard deviation (m2), skewness (m3) and
/// kurtosis (m4).  When a distribution goes long tailed, skewness and kurtosis blow up, so the
/// hope is that this type will be good for monitoring general "no bad tail" insights.
///
/// The type itelf is algebraic.  Take two readings separated by time and subtract them to get
/// perfectly recorded moments for the interval between the points.
#[derive(Clone, Default, Debug)]
pub struct Moments {
    /// The number of observations.
    pub n: u64,
    /// A value used to calculate mean.  It is mean.
    pub m1: f64,
    /// A value used to calculate variance.  It is not variance.
    pub m2: f64,
    /// A value used to calculate skewness.  It is not skewness.
    pub m3: f64,
    /// A value used to calculate kurtosis.  It is not kurtosis.
    pub m4: f64,
}

impl Moments {
    /// Create a new [Moments] with zero values.
    pub const fn new() -> Self {
        Moments {
            n: 0,
            m1: 0.0,
            m2: 0.0,
            m3: 0.0,
            m4: 0.0,
        }
    }

    /// Add `x` to the sequence of values the moments captures.
    pub fn push(&mut self, x: f64) {
        let n1: f64 = self.n as f64;
        self.n += 1;
        let n: f64 = self.n as f64;
        let delta: f64 = x - self.m1;
        let delta_n: f64 = delta / n;
        let delta_n2: f64 = delta_n * delta_n;
        let term1 = delta * delta_n * n1;
        self.m4 += term1 * delta_n2 * (n * n - 3. * n + 3.) + 6. * delta_n2 * self.m2
            - 4. * delta_n * self.m3;
        self.m3 += term1 * delta_n * (n - 2.) - 3. * delta_n * self.m2;
        self.m2 += term1;
        self.m1 += delta_n;
    }

    /// How many values does this statistical summary represent?
    pub fn n(&self) -> u64 {
        self.n
    }

    /// Average of values.
    pub fn mean(&self) -> f64 {
        self.m1
    }

    /// Square of standard deviation.
    pub fn variance(&self) -> f64 {
        self.m2 / (self.n as f64 - 1.)
    }

    /// How much the distance moves one way or another.
    pub fn skewness(&self) -> f64 {
        (self.n as f64).sqrt() * self.m3 / self.m2.powf(1.5)
    }

    /// Say something about the tail of the distribution.
    pub fn kurtosis(&self) -> f64 {
        (self.n as f64) * self.m4 / (self.m2 * self.m2)
    }

    /// Add two [Moments] together to get the equivalent of one that had the same set of values
    /// pushed.
    pub fn add(lhs: &Self, rhs: &Self) -> Self {
        let delta: f64 = rhs.m1 - lhs.m1;
        let delta2: f64 = delta * delta;
        let delta3: f64 = delta * delta2;
        let delta4: f64 = delta2 * delta2;
        let lhs_n: f64 = lhs.n as f64;
        let rhs_n: f64 = rhs.n as f64;
        let n: f64 = (lhs.n + rhs.n) as f64;
        let m1: f64 = (lhs_n * lhs.m1 + rhs_n * rhs.m1) / n;
        let m2: f64 = lhs.m2 + rhs.m2 + delta2 * lhs_n * rhs_n / n;
        let m3: f64 = lhs.m3
            + rhs.m3
            + delta3 * lhs_n * rhs_n * (lhs_n - rhs_n) / (n * n)
            + 3. * delta * (lhs_n * rhs.m2 - rhs_n * lhs.m2) / n;
        let m4: f64 = lhs.m4
            + rhs.m4
            + delta4 * lhs_n * rhs_n * (lhs_n * lhs_n - lhs_n * rhs_n + rhs_n * rhs_n)
                / (n * n * n)
            + 6. * delta2 * (lhs_n * lhs_n * rhs.m2 + rhs_n * rhs_n * lhs.m2) / (n * n)
            + 4. * delta * (lhs_n * rhs.m3 - rhs_n * lhs.m3) / n;
        Self {
            n: lhs.n + rhs.n,
            m1,
            m2,
            m3,
            m4,
        }
    }

    /// Compute lhs - rhs.
    pub fn sub(lhs: &Self, rhs: &Self) -> Self {
        let lhs_n: f64 = lhs.n as f64;
        let rhs_n: f64 = rhs.n as f64;
        let n: f64 = (rhs.n - lhs.n) as f64;
        let m1: f64 = (lhs_n * lhs.m1 - rhs_n * rhs.m1) / n;
        let delta: f64 = rhs.m1 - m1;
        let delta2: f64 = delta * delta;
        let delta3: f64 = delta * delta2;
        let delta4: f64 = delta2 * delta2;
        let m2: f64 = lhs.m2 - rhs.m2 - delta2 * rhs_n * n / lhs_n;
        let m3: f64 = lhs.m3
            - rhs.m3
            - delta3 * rhs_n * n * (n - rhs_n) / lhs_n.powf(2.)
            - 3. * delta * (n * rhs.m2 - rhs_n * m2) / lhs_n;
        let m4: f64 = lhs.m4
            - rhs.m4
            - delta4 * n * rhs_n * (n * n - n * rhs_n + rhs_n * rhs_n) / (lhs_n * lhs_n * lhs_n)
            - 6. * delta2 * (n * n * rhs.m2 + rhs_n * rhs_n * m2) / (lhs_n * lhs_n)
            - 4. * delta * (n * rhs.m3 - rhs_n * m3) / lhs_n;
        Self {
            n: lhs.n - rhs.n,
            m1,
            m2,
            m3,
            m4,
        }
    }
}

/////////////////////////////////////////////// tests //////////////////////////////////////////////

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn moments_may_be_static() {
        static MOMENTS: Moments = Moments::new();
        println!(
            "MOMENTS = {} {} {} {}",
            MOMENTS.n(),
            MOMENTS.mean(),
            MOMENTS.skewness(),
            MOMENTS.kurtosis()
        );
    }
}