src/stats/univariate/kde/kernel.rs


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

//! Kernels

use crate::stats::float::Float;

/// Kernel function
pub trait Kernel<A>: Copy + Sync
where
    A: Float,
{
    /// Apply the kernel function to the given x-value.
    fn evaluate(&self, x: A) -> A;
}

/// Gaussian kernel
#[derive(Clone, Copy)]
pub struct Gaussian;

impl<A> Kernel<A> for Gaussian
where
    A: Float,
{
    fn evaluate(&self, x: A) -> A {
        use std::f32::consts::PI;

        (x.powi(2).exp() * A::cast(2. * PI)).sqrt().recip()
    }
}

#[cfg(test)]
macro_rules! test {
    ($ty:ident) => {
        mod $ty {
            mod gaussian {
                use approx::relative_eq;
                use quickcheck::quickcheck;
                use quickcheck::TestResult;

                use crate::stats::univariate::kde::kernel::{Gaussian, Kernel};

                quickcheck! {
                    fn symmetric(x: $ty) -> bool {
                        relative_eq!(Gaussian.evaluate(-x), Gaussian.evaluate(x))
                    }
                }

                // Any [a b] integral should be in the range [0 1]
                quickcheck! {
                    fn integral(a: $ty, b: $ty) -> TestResult {
                        const DX: $ty = 1e-3;

                        if a > b {
                            TestResult::discard()
                        } else {
                            let mut acc = 0.;
                            let mut x = a;
                            let mut y = Gaussian.evaluate(a);

                            while x < b {
                                acc += DX * y / 2.;

                                x += DX;
                                y = Gaussian.evaluate(x);

                                acc += DX * y / 2.;
                            }

                            TestResult::from_bool(
                                (acc > 0. || relative_eq!(acc, 0.)) &&
                                (acc < 1. || relative_eq!(acc, 1.)))
                        }
                    }
                }
            }
        }
    };
}

#[cfg(test)]
mod test {
    test!(f32);
    test!(f64);
}