statistical/
univariate_.rs

1// Copyright (c) 2015 Jeff Belgum
2//
3// Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated
4// documentation files (the "Software"), to deal in the Software without restriction, including without
5// limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of
6// the Software, and to permit persons to whom the Software is furnished to do so, subject to the following
7// conditions:
8//
9// The above copyright notice and this permission notice shall be included in all copies or substantial portions
10// of the Software.
11//
12// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED
13// TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT
14// SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
15// OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
16// DEALINGS IN THE SOFTWARE.
17
18extern crate rand;
19extern crate num;
20
21use std::collections::HashMap;
22use std::hash::Hash;
23
24use num::{Float,
25          One,
26          PrimInt,
27          Zero};
28
29use super::stats_ as stats;
30
31pub fn harmonic_mean<T>(v: &[T]) -> T
32    where T: Float
33{
34    let invert = |x: &T| T::one() / *x;
35    let sum_of_inverted = v.iter().map(invert).fold(T::zero(), |acc, elem| acc + elem);
36    num::cast::<usize, T>(v.len()).unwrap() / sum_of_inverted
37}
38
39pub fn geometric_mean<T>(v: &[T]) -> T
40    where T: Float
41{
42    let product = v.iter().fold(T::one(), |acc, elem| acc * *elem);
43    let one_over_len = T::one() / num::cast(v.len()).unwrap();
44    product.powf(one_over_len)
45}
46
47pub fn quadratic_mean<T>(v: &[T]) -> T
48    where T: Float
49{
50    let square = |x: &T| (*x).powi(2);
51    let sum_of_squared = v.iter().map(square).fold(T::zero(), |acc, elem| acc + elem);
52    (sum_of_squared / num::cast(v.len()).unwrap()).sqrt()
53}
54
55pub fn mode<T>(v: &[T]) -> Option<T>
56    where T: Hash + Copy + Eq
57{
58    match v.len() {
59        0 => None,
60        1 => Some(v[0]),
61        _ => {
62            let mut counter = HashMap::new();
63            for x in v.iter() {
64                let count = counter.entry(x).or_insert(0);
65                *count += 1;
66            }
67            let mut max = -1;
68            let mut mode = None;
69
70            for (val, count) in counter.iter() {
71                if *count > max {
72                    max = *count;
73                    mode = Some(**val);
74                }
75            }
76            mode
77        }
78    }
79}
80
81pub fn average_deviation<T>(v: &[T], mean: Option<T>) -> T
82    where T: Float
83{
84    let mean = mean.unwrap_or_else(|| stats::mean(v));
85    let dev = v.iter().map(|&x| (x-mean).abs()).fold(T::zero(), |acc, elem| acc + elem);
86    dev / num::cast(v.len()).unwrap()
87
88}
89
90pub fn pearson_skewness<T>(mean: T, mode: T, stdev: T) -> T
91    where T: Float
92{
93    (mean - mode) / stdev
94}
95
96pub fn skewness<T>(v: &[T], mean: Option<T>, pstdev: Option<T>) -> T
97    where T: Float
98{
99    let m = stats::std_moment(v, stats::Degree::Three, mean, pstdev);
100    let n = num::cast(v.len()).unwrap();
101    let skew = m / n;
102    let k = ( n * ( n - T::one())).sqrt()/( n - num::cast(2).unwrap());
103    skew * k
104}
105
106pub fn pskewness<T>(v: &[T], mean: Option<T>, pstdev: Option<T>) -> T
107    where T: Float
108{
109    let m = stats::std_moment(v, stats::Degree::Three, mean, pstdev);
110    m / num::cast(v.len()).unwrap()
111}
112
113pub fn kurtosis<T>(v: &[T], mean: Option<T>, pstdev: Option<T>) -> T
114    where T: Float
115{
116    let two = num::cast::<f32, T>(2.0).unwrap();
117    let three = num::cast::<f32, T>(3.0).unwrap();
118
119    let m = stats::std_moment(v, stats::Degree::Four, mean, pstdev);
120    let n = num::cast(v.len()).unwrap();
121    let q = (n - T::one())/((n-two)*(n-three));
122    let gamma2 = m / n;
123    let kurt = q * (( ( n + T::one() ) * gamma2) - ( (n-T::one()) * three ));
124    kurt
125}
126
127pub fn pkurtosis<T>(v: &[T], mean: Option<T>, pstdev: Option<T>) -> T
128    where T: Float
129{
130    let m = stats::std_moment(v, stats::Degree::Four, mean, pstdev);
131    m / num::cast(v.len()).unwrap() - num::cast(3).unwrap()
132}
133
134pub fn standard_error_mean<T>(stdev: T, sample_size: T, population_size: Option<T>) -> T
135    where T: Float
136{
137    let mut err = stdev / sample_size.sqrt();
138    if let Some(p) = population_size {
139        err = err * ((p - sample_size) / (p - T::one())).sqrt()
140    }
141    err
142
143}
144
145pub fn standard_error_skewness<T, U>(sample_size: T) -> U
146    where T: PrimInt, U: Float
147{
148    (num::cast::<f32,U>(6.0).unwrap() / num::cast(sample_size).unwrap()).sqrt()
149}
150
151pub fn standard_error_kurtosis<T, U>(sample_size: T) -> U
152    where T: PrimInt, U: Float
153{
154    (num::cast::<f32,U>(24.0).unwrap() / num::cast(sample_size).unwrap()).sqrt()
155}
156
157#[cfg(test)]
158mod test {
159    use super::*;
160
161    #[test]
162    fn test_harmonic_mean() {
163        let vec = vec![0.25, 0.5, 1.0, 1.0];
164        assert_eq!(harmonic_mean(&vec), 0.5);
165        let vec = vec![0.5, 0.5, 0.5];
166        assert_eq!(harmonic_mean(&vec), 0.5);
167        let vec = vec![1.0,2.0,4.0];
168        assert_eq!(harmonic_mean(&vec), 12.0/7.0);
169    }
170    #[test]
171    fn test_geometric_mean() {
172        let vec = vec![1.0, 2.0, 6.125, 12.25];
173        assert_eq!(geometric_mean(&vec), 3.5);
174    }
175    #[test]
176    fn test_quadratic_mean() {
177        let vec = vec![-3.0, -2.0, 0.0, 2.0, 3.0];
178        assert_eq!(quadratic_mean(&vec), 2.280350850198276);
179    }
180    #[test]
181    fn test_mode() {
182        let vec = vec![2,4,3,5,4,6,1,1,6,4,0,0];
183        assert_eq!(mode(&vec), Some(4));
184        let vec = vec![1];
185        assert_eq!(mode(&vec), Some(1));
186    }
187    #[test]
188    fn test_average_deviation() {
189        let vec = vec![2.0, 2.25, 2.5, 2.5, 3.25];
190        assert_eq!(average_deviation(&vec, None), 0.3);
191        assert_eq!(average_deviation(&vec, Some(2.75)), 0.45);
192    }
193    #[test]
194    fn test_pearson_skewness() {
195        assert_eq!(pearson_skewness(2.5, 2.25, 2.5), 0.1);
196        assert_eq!(pearson_skewness(2.5, 5.75, 0.5), -6.5);
197    }
198    #[test]
199    fn test_skewness() {
200        let vec = vec![1.25, 1.5, 1.5, 1.75, 1.75, 2.5, 2.75, 4.5];
201        assert_eq!(skewness(&vec, None, None), 1.7146101353987853);
202        let vec = vec![1.25, 1.5, 1.5, 1.75, 1.75, 2.5, 2.75, 4.5];
203        assert_eq!(skewness(&vec, Some(2.25), Some(1.0)), 1.4713288161532945);
204    }
205    #[test]
206    fn test_pskewness() {
207        let vec = vec![1.25, 1.5, 1.5, 1.75, 1.75, 2.5, 2.75, 4.5];
208        assert_eq!(pskewness(&vec, None, None), 1.3747465025469285);
209    }
210    #[test]
211    fn test_kurtosis() {
212        let vec = vec![1.25, 1.5, 1.5, 1.75, 1.75, 2.5, 2.75, 4.5];
213        assert_eq!(kurtosis(&vec, None, None), 3.036788927335642);
214        let vec = vec![1.25, 1.5, 1.5, 1.75, 1.75, 2.5, 2.75, 4.5];
215        assert_eq!(kurtosis(&vec, Some(2.25), Some(1.0)), 2.3064453125);
216    }
217    #[test]
218    fn test_pkurtosis() {
219        let vec = vec![1.25, 1.5, 1.5, 1.75, 1.75, 2.5, 2.75, 4.5];
220        assert_eq!(pkurtosis(&vec, None, None), 0.7794232987312579);
221    }
222    #[test]
223    fn test_standard_error_mean() {
224        assert_eq!(standard_error_mean(2.0, 16.0, None), 0.5);
225    }
226    #[test]
227    fn test_standard_error_skewness() {
228        assert_eq!(standard_error_skewness::<i32, f32>(15), 0.63245553203);
229    }
230    #[test]
231    fn test_standard_error_kurtosis() {
232        assert_eq!(standard_error_kurtosis::<i32, f32>(15), 1.2649110640);
233    }
234}
statistical/univariate_.rs

statistical/
univariate_.rs
