criterion/stats/univariate/kde/
kernel.rs

1//! Kernels
2
3use crate::stats::float::Float;
4
5/// Kernel function
6pub trait Kernel<A>: Copy + Sync
7where
8    A: Float,
9{
10    /// Apply the kernel function to the given x-value.
11    fn evaluate(&self, x: A) -> A;
12}
13
14/// Gaussian kernel
15#[derive(Clone, Copy)]
16pub struct Gaussian;
17
18impl<A> Kernel<A> for Gaussian
19where
20    A: Float,
21{
22    fn evaluate(&self, x: A) -> A {
23        use std::f32::consts::PI;
24
25        (x.powi(2).exp() * A::cast(2. * PI)).sqrt().recip()
26    }
27}
28
29#[cfg(test)]
30macro_rules! test {
31    ($ty:ident) => {
32        mod $ty {
33            mod gaussian {
34                use approx::relative_eq;
35                use quickcheck::quickcheck;
36                use quickcheck::TestResult;
37
38                use crate::stats::univariate::kde::kernel::{Gaussian, Kernel};
39
40                quickcheck! {
41                    fn symmetric(x: $ty) -> bool {
42                        x.is_nan() || relative_eq!(Gaussian.evaluate(-x), Gaussian.evaluate(x))
43                    }
44                }
45
46                // Any [a b] integral should be in the range [0 1]
47                quickcheck! {
48                    fn integral(a: $ty, b: $ty) -> TestResult {
49                        let a = a.sin().abs(); // map the value to [0 1]
50                        let b = b.sin().abs(); // map the value to [0 1]
51                        const DX: $ty = 1e-3;
52
53                        if a > b {
54                            TestResult::discard()
55                        } else {
56                            let mut acc = 0.;
57                            let mut x = a;
58                            let mut y = Gaussian.evaluate(a);
59
60                            while x < b {
61                                acc += DX * y / 2.;
62
63                                x += DX;
64                                y = Gaussian.evaluate(x);
65
66                                acc += DX * y / 2.;
67                            }
68
69                            TestResult::from_bool(
70                                (acc > 0. || relative_eq!(acc, 0.)) &&
71                                (acc < 1. || relative_eq!(acc, 1.)))
72                        }
73                    }
74                }
75            }
76        }
77    };
78}
79
80#[cfg(test)]
81mod test {
82    test!(f32);
83    test!(f64);
84}