statest 0.2.2

Rust crate for statistical test.
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

//! Kolmogorov–Smirnov test.

use ndarray::*;
use num_traits::Float;
use statrs::distribution::Univariate;

use crate::cast::*;

/// Max iteration when calculating Q_KS.
const MAXJ:usize = 100;
/// Small value 1.
const EPS1:f64 = 1.0e-3;
/// Small value 2.
const EPS2:f64 = 1.0e-8;

/// Struct for Kolmogorov–Smirnov test.
pub struct KSTest {
    pub n: usize,
    x: Vec<f64>,
}

impl KSTest {
    pub fn new(x: &[f64]) -> Self {
        let n = x.len();
        let mut x = x.to_vec();
        x.sort_by(|&a, b| a.partial_cmp(b).unwrap());
        Self { n, x }
    }

    /// Kolmogorov–Smirnov test which returns probability `prob` and test statistic D `d`. 
    pub fn ks1<T: Univariate<f64, f64>>(&self, dist: &T) -> (f64, f64) {
        let mut d = 0.0;
        let mut f_old = 0.0;

        for (i, &xi) in self.x.iter().enumerate() {
            let f_new = (i + 1) as f64 / self.n as f64;
            let f_dist = dist.cdf(xi);
            let (df_old, df_new) = ((f_old - f_dist).abs(), (f_new - f_dist).abs());
            let dt = if df_old > df_new {
                df_old
            } else {
                df_new
            };
            if dt > d {
                d = dt;
            }
            f_old = f_new;
        }
        let n_sqrt = (self.n as f64).sqrt();
        let prob = q_ks((n_sqrt + 0.12 + 0.11 / n_sqrt) * d);
        (prob, d)
    }
}

/// Trait for Kolmogorov–Smirnov test.
pub trait KSVec<T: Float> {
    /// Kolmogorov–Smirnov test. If `p` <= `prob`, then returns `true`.
    fn ks1<S: Univariate<f64, f64>>(&self, dist: &S, p: T) -> bool;
}

impl KSVec<f64> for Vec<f64> {
    fn ks1<S: Univariate<f64, f64>>(&self, dist: &S, p: f64) -> bool {
        let ks = KSTest::new(&self);
        let (prob, _) = ks.ks1(dist);
        prob >= p
    }
}

impl KSVec<f32> for Vec<f32> {
    fn ks1<S: Univariate<f64, f64>>(&self, dist: &S, p: f32) -> bool {
        self.to_vec_f64().ks1(dist, p as f64)
    }
}

impl<U: Data<Elem = f64>> KSVec<f64> for ArrayBase<U, Ix1> {
    fn ks1<S: Univariate<f64, f64>>(&self, dist: &S, p: f64) -> bool {
        self.to_vec().ks1(dist, p)
    }
}

impl<U: Data<Elem = f32>> KSVec<f32> for ArrayBase<U, Ix1> {
    fn ks1<S: Univariate<f64, f64>>(&self, dist: &S, p: f32) -> bool {
        self.to_vec().ks1(dist, p)
    }
}

/// Calculate test statistic Q_KS.
fn q_ks(lambda: f64) -> f64 {
    let mut sum = 0.0;
    let mut term_old = 0.0;
    let mut fac = 2.0;
    let m2lambda2 = -2.0 * lambda * lambda;
    for i in 1..=MAXJ {
        let j = i as f64;
        let term = fac * (m2lambda2 * j * j).exp();
        sum += term;
        if term.abs() <= EPS1 * term_old || term.abs() <= EPS2 * sum {
            return sum;
        }
        fac = -fac;
        term_old = term.abs();
    }
    1.0
}