average/moments/
kurtosis.rs

1/// Estimate the arithmetic mean, the variance, the skewness and the kurtosis of
2/// a sequence of numbers ("population").
3///
4/// This can be used to estimate the standard error of the mean.
5#[derive(Debug, Clone)]
6#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
7pub struct Kurtosis {
8    /// Estimator of mean, variance and skewness.
9    avg: Skewness,
10    /// Intermediate sum of terms to the fourth for calculating the skewness.
11    sum_4: f64,
12}
13
14impl Kurtosis {
15    /// Create a new kurtosis estimator.
16    #[inline]
17    pub fn new() -> Kurtosis {
18        Kurtosis {
19            avg: Skewness::new(),
20            sum_4: 0.,
21        }
22    }
23
24    /// Increment the sample size.
25    ///
26    /// This does not update anything else.
27    #[inline]
28    fn increment(&mut self) {
29        self.avg.increment();
30    }
31
32    /// Add an observation given an already calculated difference from the mean
33    /// divided by the number of samples, assuming the inner count of the sample
34    /// size was already updated.
35    ///
36    /// This is useful for avoiding unnecessary divisions in the inner loop.
37    #[inline]
38    fn add_inner(&mut self, delta: f64, delta_n: f64) {
39        // This algorithm was suggested by Terriberry.
40        //
41        // See https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance.
42        let n = self.len().to_f64().unwrap();
43        let term = delta * delta_n * (n - 1.);
44        let delta_n_sq = delta_n*delta_n;
45        self.sum_4 += term * delta_n_sq * (n*n - 3.*n + 3.)
46            + 6. * delta_n_sq * self.avg.avg.sum_2
47            - 4. * delta_n * self.avg.sum_3;
48        self.avg.add_inner(delta, delta_n);
49    }
50
51    /// Determine whether the sample is empty.
52    #[inline]
53    pub fn is_empty(&self) -> bool {
54        self.avg.is_empty()
55    }
56
57    /// Estimate the mean of the population.
58    ///
59    /// Returns NaN for an empty sample.
60    #[inline]
61    pub fn mean(&self) -> f64 {
62        self.avg.mean()
63    }
64
65    /// Return the sample size.
66    #[inline]
67    pub fn len(&self) -> u64 {
68        self.avg.len()
69    }
70
71    /// Calculate the sample variance.
72    ///
73    /// This is an unbiased estimator of the variance of the population.
74    /// 
75    /// Returns NaN for samples of size 1 or less.
76    #[inline]
77    pub fn sample_variance(&self) -> f64 {
78        self.avg.sample_variance()
79    }
80
81    /// Calculate the population variance of the sample.
82    ///
83    /// This is a biased estimator of the variance of the population.
84    /// 
85    /// Returns NaN for an empty sample.
86    #[inline]
87    pub fn population_variance(&self) -> f64 {
88        self.avg.population_variance()
89    }
90
91    /// Estimate the standard error of the mean of the population.
92    #[inline]
93    pub fn error_mean(&self) -> f64 {
94        self.avg.error_mean()
95    }
96
97    /// Estimate the skewness of the population.
98    /// 
99    /// Returns NaN for an empty sample.
100    #[inline]
101    pub fn skewness(&self) -> f64 {
102        self.avg.skewness()
103    }
104
105    /// Estimate the excess kurtosis of the population.
106    /// 
107    /// Returns NaN for an empty sample.
108    #[inline]
109    pub fn kurtosis(&self) -> f64 {
110        if self.is_empty() {
111            return f64::NAN;
112        }
113        if self.sum_4 == 0. {
114            return 0.;
115        }
116        let n = self.len().to_f64().unwrap();
117        debug_assert_ne!(self.avg.avg.sum_2, 0.);
118        n * self.sum_4 / (self.avg.avg.sum_2 * self.avg.avg.sum_2) - 3.
119    }
120
121}
122
123impl core::default::Default for Kurtosis {
124    fn default() -> Kurtosis {
125        Kurtosis::new()
126    }
127}
128
129impl Estimate for Kurtosis {
130    #[inline]
131    fn add(&mut self, x: f64) {
132        let delta = x - self.avg.avg.avg.avg;
133        self.increment();
134        let n = self.len().to_f64().unwrap();
135        self.add_inner(delta, delta/n);
136    }
137
138    #[inline]
139    fn estimate(&self) -> f64 {
140        self.kurtosis()
141    }
142}
143
144impl Merge for Kurtosis {
145    #[inline]
146    fn merge(&mut self, other: &Kurtosis) {
147        if other.is_empty() {
148            return;
149        }
150        if self.is_empty() {
151            *self = other.clone();
152            return;
153        }
154        let len_self = self.len().to_f64().unwrap();
155        let len_other = other.len().to_f64().unwrap();
156        let len_total = len_self + len_other;
157        let delta = other.mean() - self.mean();
158        let delta_n = delta / len_total;
159        let delta_n_sq = delta_n * delta_n;
160        self.sum_4 += other.sum_4
161            + delta * delta_n*delta_n_sq * len_self*len_other
162              * (len_self*len_self - len_self*len_other + len_other*len_other)
163            + 6.*delta_n_sq * (len_self*len_self * other.avg.avg.sum_2 + len_other*len_other * self.avg.avg.sum_2)
164            + 4.*delta_n * (len_self * other.avg.sum_3 - len_other * self.avg.sum_3);
165        self.avg.merge(&other.avg);
166    }
167}
168
169impl_from_iterator!(Kurtosis);
170impl_from_par_iterator!(Kurtosis);
171impl_extend!(Kurtosis);
average/moments/kurtosis.rs

average/moments/
kurtosis.rs
