Add more Kolmogorov–Smirnov test

rand_distr/tests/cdf.rs

@@ -11,8 +11,7 @@ use core::f64;
 use num_traits::AsPrimitive;
 use rand::SeedableRng;
 use rand_distr::{Distribution, Normal};
-use special::Beta;
-use special::Primitive;
+use special::{Beta, Gamma, Primitive};
 
 // [1] Nonparametric Goodness-of-Fit Tests for Discrete Null Distributions
 //     by Taylor B. Arnold and John W. Emerson
@@ -159,6 +158,186 @@ fn normal() {
     }
 }
 
+#[test]
+fn cauchy() {
+    fn cdf(x: f64, median: f64, scale: f64) -> f64 {
+        (1.0 / f64::consts::PI) * ((x - median) / scale).atan() + 0.5
+    }
+
+    let parameters = [
+        (0.0, 1.0),
+        (0.0, 0.1),
+        (1.0, 10.0),
+        (1.0, 100.0),
+        (-1.0, 0.00001),
+        (-1.0, 0.0000001),
+    ];
+
+    for (seed, (median, scale)) in parameters.into_iter().enumerate() {
+        let dist = rand_distr::Cauchy::new(median, scale).unwrap();
+        test_continuous(seed as u64, dist, |x| cdf(x, median, scale));
+    }
+}
+
+#[test]
+fn uniform() {
+    fn cdf(x: f64, a: f64, b: f64) -> f64 {
+        if x < a {
+            0.0
+        } else if x < b {
+            (x - a) / (b - a)
+        } else {
+            1.0
+        }
+    }
+
+    let parameters = [(0.0, 1.0), (-1.0, 1.0), (0.0, 100.0), (-100.0, 100.0)];
+
+    for (seed, (a, b)) in parameters.into_iter().enumerate() {
+        let dist = rand_distr::Uniform::new(a, b).unwrap();
+        test_continuous(seed as u64, dist, |x| cdf(x, a, b));
+    }
+}
+
+#[test]
+fn log_normal() {
+    fn cdf(x: f64, mean: f64, std_dev: f64) -> f64 {
+        if x <= 0.0 {
+            0.0
+        } else if x.is_infinite() {
+            1.0
+        } else {
+            0.5 * ((mean - x.ln()) / (std_dev * f64::consts::SQRT_2)).erfc()
+        }
+    }
+
+    let parameters = [
+        (0.0, 1.0),
+        (0.0, 0.1),
+        (0.5, 0.7),
+        (1.0, 10.0),
+        (1.0, 100.0),
+    ];
+
+    println!("{}", cdf(5.0, 2.0, 1.5));
+
+    for (seed, (mean, std_dev)) in parameters.into_iter().enumerate() {
+        let dist = rand_distr::LogNormal::new(mean, std_dev).unwrap();
+        test_continuous(seed as u64, dist, |x| cdf(x, mean, std_dev));
+    }
+}
+
+#[test]
+fn exp() {
+    fn cdf(x: f64, lambda: f64) -> f64 {
+        1.0 - (-lambda * x).exp()
+    }
+
+    let parameters = [0.5, 1.0, 7.5, 32.0, 100.0];
+
+    for (seed, lambda) in parameters.into_iter().enumerate() {
+        let dist = rand_distr::Exp::new(lambda).unwrap();
+        test_continuous(seed as u64, dist, |x| cdf(x, lambda));
+    }
+}
+
+#[test]
+fn weibull() {
+    fn cdf(x: f64, lambda: f64, k: f64) -> f64 {
+        if x < 0.0 {
+            return 0.0;
+        }
+
+        1.0 - (-(x / lambda).powf(k)).exp()
+    }
+
+    let parameters = [
+        (0.5, 1.0),
+        (1.0, 1.0),
+        (10.0, 0.1),
+        (0.1, 10.0),
+        (15.0, 20.0),
+        // (1000.0, 0.001), // Fail case
+    ];
+
+    for (seed, (lambda, k)) in parameters.into_iter().enumerate() {
+        let dist = rand_distr::Weibull::new(lambda, k).unwrap();
+        test_continuous(seed as u64, dist, |x| cdf(x, lambda, k));
+    }
+}
+
+#[test]
+fn gamma() {
+    fn cdf(x: f64, shape: f64, scale: f64) -> f64 {
+        if x < 0.0 {
+            return 0.0;
+        }
+
+        (x / scale).inc_gamma(shape)
+    }
+
+    let parameters = [
+        (0.5, 2.0),
+        (1.0, 1.0),
+        (10.0, 0.1),
+        (100.0, 0.0001),
+        (0.9999, 2.0),
+    ];
+
+    for (seed, (shape, scale)) in parameters.into_iter().enumerate() {
+        let dist = rand_distr::Gamma::new(shape, scale).unwrap();
+        test_continuous(seed as u64, dist, |x| cdf(x, shape, scale));
+    }
+}
+
+#[test]
+fn chi_squared() {
+    fn cdf(x: f64, k: f64) -> f64 {
+        if x < 0.0 {
+            return 0.0;
+        }
+
+        (x / 2.0).inc_gamma(k / 2.0)
+    }
+
+    let parameters = [
+        // 0.01, // Fail case
+        0.1, 1.0, 2.0, 10.0, 100.0, 1000.0,
+    ];
+
+    for (seed, k) in parameters.into_iter().enumerate() {
+        let dist = rand_distr::ChiSquared::new(k).unwrap();
+        test_continuous(seed as u64, dist, |x| cdf(x, k));
+    }
+}
+
+#[test]
+fn beta() {
+    fn cdf(x: f64, alpha: f64, beta: f64) -> f64 {
+        if x < 0.0 {
+            return 0.0;
+        }
+        if x > 1.0 {
+            return 1.0;
+        }
+        let ln_beta_ab = alpha.ln_beta(beta);
+        x.inc_beta(alpha, beta, ln_beta_ab)
+    }
+
+    let parameters = [
+        (0.5, 0.5),
+        (2.0, 3.5),
+        (10.0, 1.0),
+        (100.0, 50.0),
+        // (10.0, 0.1), // Fail case
+    ];
+
+    for (seed, (alpha, beta)) in parameters.into_iter().enumerate() {
+        let dist = rand_distr::Beta::new(alpha, beta).unwrap();
+        test_continuous(seed as u64, dist, |x| cdf(x, alpha, beta));
+    }
+}
+
 fn binomial_cdf(k: i64, p: f64, n: u64) -> f64 {
     if k < 0 {
         return 0.0;
@@ -195,3 +374,81 @@ fn binomial() {
         });
     }
 }
+
+#[test]
+fn geometric() {
+    fn cdf(k: i64, p: f64) -> f64 {
+        if k < 0 {
+            0.0
+        } else {
+            1.0 - (1.0 - p).powi(1 + k as i32)
+        }
+    }
+
+    let parameters = [0.3, 0.5, 0.7, 0.0000001, 0.9999];
+
+    for (seed, p) in parameters.into_iter().enumerate() {
+        let dist = rand_distr::Geometric::new(p).unwrap();
+        test_discrete(seed as u64, dist, |k| cdf(k, p));
+    }
+}
+
+#[test]
+fn hypergeometric() {
+    fn cdf(x: i64, n: u64, k: u64, n_: u64) -> f64 {
+        let min = if n_ + k > n { n_ + k - n } else { 0 };
+        let max = k.min(n_);
+        if x < min as i64 {
+            return 0.0;
+        } else if x >= max as i64 {
+            return 1.0;
+        }
+
+        (min..x as u64 + 1).fold(0.0, |acc, k_| {
+            acc + (ln_binomial(k, k_) + ln_binomial(n - k, n_ - k_) - ln_binomial(n, n_)).exp()
+        })
+    }
+
+    let parameters = [
+        (15, 13, 10),
+        (25, 15, 5),
+        (60, 10, 7),
+        (70, 20, 50),
+        (100, 50, 10),
+        // (100, 50, 49), // Fail case
+    ];
+
+    for (seed, (n, k, n_)) in parameters.into_iter().enumerate() {
+        let dist = rand_distr::Hypergeometric::new(n, k, n_).unwrap();
+        test_discrete(seed as u64, dist, |x| cdf(x, n, k, n_));
+    }
+}
+
+#[test]
+fn poisson() {
+    fn cdf(k: i64, lambda: f64) -> f64 {
+        if k < 0 || lambda <= 0.0 {
+            return 0.0;
+        }
+
+        1.0 - lambda.inc_gamma(k as f64 + 1.0)
+    }
+
+    let parameters = [
+        0.1_f32, 1.0, 7.5,
+        // 1.844E+19,  // fail case
+    ];
+
+    for (seed, lambda) in parameters.into_iter().enumerate() {
+        let dist = rand_distr::Poisson::new(lambda).unwrap();
+        test_discrete(seed as u64, dist, |k| cdf(k, lambda as f64));
+    }
+}
+
+fn ln_factorial(n: u64) -> f64 {
+    (n as f64 + 1.0).lgamma().0
+}
+
+fn ln_binomial(n: u64, k: u64) -> f64 {
+    ln_factorial(n) - ln_factorial(k) - ln_factorial(n - k)
+}