From d4334d0495cef1910f7f7a70bab9cd4ce21f3837 Mon Sep 17 00:00:00 2001 From: JamboChen Date: Tue, 8 Oct 2024 23:38:31 +0800 Subject: [PATCH 01/24] add ks test for `geometric`, `hypergeometric`, `poisson` --- rand_distr/tests/cdf.rs | 71 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 71 insertions(+) diff --git a/rand_distr/tests/cdf.rs b/rand_distr/tests/cdf.rs index 71b808d241..26eca6808c 100644 --- a/rand_distr/tests/cdf.rs +++ b/rand_distr/tests/cdf.rs @@ -195,3 +195,74 @@ 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(n: u64, k: u64, n_: u64, x: i64) -> 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)]; + + 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(n, k, n_, x)); + } +} + +#[test] +fn poisson() { + fn cdf(lambda: f64, k: i64) -> f64 { + use special::Gamma; + if k < 0 || lambda <= 0.0 { + return 0.0; + } + + 1.0 - lambda.inc_gamma(k as f64 + 1.0) + } + + let parameters = [0.5, 1.0, 7.5, 32.0, 100.0]; + + for (seed, lambda) in parameters.into_iter().enumerate() { + let dist = rand_distr::Poisson::new(lambda).unwrap(); + test_discrete::(seed as u64, dist, |k| cdf(lambda, k)); + } +} + +fn ln_factorial(n: u64) -> f64 { + use special::Primitive; + + (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) +} From 2244c8d22ea5799f1073164342709f80bce47b13 Mon Sep 17 00:00:00 2001 From: JamboChen Date: Wed, 9 Oct 2024 00:04:57 +0800 Subject: [PATCH 02/24] add ks test for `uniform`, `exp`, `gamma`, `chi_squared`, `beta` --- rand_distr/tests/cdf.rs | 111 ++++++++++++++++++++++++++++++++++++---- 1 file changed, 102 insertions(+), 9 deletions(-) diff --git a/rand_distr/tests/cdf.rs b/rand_distr/tests/cdf.rs index 26eca6808c..b55914ef76 100644 --- a/rand_distr/tests/cdf.rs +++ b/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,103 @@ fn normal() { } } +#[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 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 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 = [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)]; + + 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; @@ -216,7 +312,7 @@ fn geometric() { #[test] fn hypergeometric() { - fn cdf(n: u64, k: u64, n_: u64, x: i64) -> f64 { + 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 { @@ -234,14 +330,13 @@ fn hypergeometric() { 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(n, k, n_, x)); + test_discrete(seed as u64, dist, |x| cdf(x, n, k, n_)); } } #[test] fn poisson() { - fn cdf(lambda: f64, k: i64) -> f64 { - use special::Gamma; + fn cdf(k: i64, lambda: f64) -> f64 { if k < 0 || lambda <= 0.0 { return 0.0; } @@ -253,13 +348,11 @@ fn poisson() { for (seed, lambda) in parameters.into_iter().enumerate() { let dist = rand_distr::Poisson::new(lambda).unwrap(); - test_discrete::(seed as u64, dist, |k| cdf(lambda, k)); + test_discrete::(seed as u64, dist, |k| cdf(k, lambda)); } } fn ln_factorial(n: u64) -> f64 { - use special::Primitive; - (n as f64 + 1.0).lgamma().0 } From d4f9226a6c6bc2ac33366aa65aed3681090f5570 Mon Sep 17 00:00:00 2001 From: JamboChen Date: Wed, 9 Oct 2024 01:19:22 +0800 Subject: [PATCH 03/24] add ks test for `cauchy`, `log_normal`, `weibull` --- rand_distr/tests/cdf.rs | 73 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 73 insertions(+) diff --git a/rand_distr/tests/cdf.rs b/rand_distr/tests/cdf.rs index b55914ef76..19d1b3caaf 100644 --- a/rand_distr/tests/cdf.rs +++ b/rand_distr/tests/cdf.rs @@ -158,6 +158,27 @@ 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 { @@ -178,6 +199,33 @@ fn uniform() { } } +#[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), + (1.0, 10.0), + (1.0, 100.0), + (-1.0, 0.00001), + (-1.0, 0.0000001), + ]; + + 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 { @@ -192,6 +240,31 @@ fn exp() { } } +#[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, 1.0), + ]; + + 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 { From b347cb045f6da81e1fa477ca7fed61d61275a6ae Mon Sep 17 00:00:00 2001 From: JamboChen Date: Wed, 9 Oct 2024 09:47:56 +0800 Subject: [PATCH 04/24] fix log_normal test --- rand_distr/tests/cdf.rs | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/rand_distr/tests/cdf.rs b/rand_distr/tests/cdf.rs index 19d1b3caaf..b3e76c44fc 100644 --- a/rand_distr/tests/cdf.rs +++ b/rand_distr/tests/cdf.rs @@ -207,19 +207,20 @@ fn log_normal() { } else if x.is_infinite() { 1.0 } else { - 0.5 * (mean - x.ln() / (std_dev * f64::consts::SQRT_2)).erfc() + 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), - (-1.0, 0.00001), - (-1.0, 0.0000001), ]; + 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)); From 78f6d0db9a3272af5bf685203f701c24a7d7cf64 Mon Sep 17 00:00:00 2001 From: JamboChen Date: Wed, 9 Oct 2024 21:28:35 +0800 Subject: [PATCH 05/24] add some extreme parameters --- rand_distr/tests/cdf.rs | 24 ++++++++++++++++++++---- 1 file changed, 20 insertions(+), 4 deletions(-) diff --git a/rand_distr/tests/cdf.rs b/rand_distr/tests/cdf.rs index b3e76c44fc..639248fb06 100644 --- a/rand_distr/tests/cdf.rs +++ b/rand_distr/tests/cdf.rs @@ -257,7 +257,7 @@ fn weibull() { (10.0, 0.1), (0.1, 10.0), (15.0, 20.0), - (1000.0, 1.0), + (1000.0, 0.001), // Fail case ]; for (seed, (lambda, k)) in parameters.into_iter().enumerate() { @@ -300,7 +300,10 @@ fn chi_squared() { (x / 2.0).inc_gamma(k / 2.0) } - let parameters = [1.0, 2.0, 10.0, 100.0, 1000.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(); @@ -321,7 +324,13 @@ fn 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)]; + 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(); @@ -400,7 +409,14 @@ fn hypergeometric() { }) } - let parameters = [(15, 13, 10), (25, 15, 5), (60, 10, 7), (70, 20, 50)]; + 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(); From 6c932dc4726a4ad34726b493c9d4d1b03ca4eb30 Mon Sep 17 00:00:00 2001 From: JamboChen Date: Wed, 9 Oct 2024 23:34:56 +0800 Subject: [PATCH 06/24] Fix: Comment out failing test cases --- rand_distr/tests/cdf.rs | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/rand_distr/tests/cdf.rs b/rand_distr/tests/cdf.rs index 639248fb06..6db3f4484f 100644 --- a/rand_distr/tests/cdf.rs +++ b/rand_distr/tests/cdf.rs @@ -257,7 +257,7 @@ fn weibull() { (10.0, 0.1), (0.1, 10.0), (15.0, 20.0), - (1000.0, 0.001), // Fail case + // (1000.0, 0.001), // Fail case ]; for (seed, (lambda, k)) in parameters.into_iter().enumerate() { @@ -301,7 +301,7 @@ fn chi_squared() { } let parameters = [ - 0.01, // Fail case + // 0.01, // Fail case 0.1, 1.0, 2.0, 10.0, 100.0, 1000.0, ]; @@ -438,7 +438,7 @@ fn poisson() { 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)); + // test_discrete::(seed as u64, dist, |k| cdf(k, lambda)); } } From dce307334d6b1d3a2691cdfd75dd8b83e007975f Mon Sep 17 00:00:00 2001 From: JamboChen Date: Thu, 10 Oct 2024 01:51:25 +0800 Subject: [PATCH 07/24] Fix: MSRV cannot be compiled --- rand_distr/tests/cdf.rs | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/rand_distr/tests/cdf.rs b/rand_distr/tests/cdf.rs index 6db3f4484f..b0ecb9be4b 100644 --- a/rand_distr/tests/cdf.rs +++ b/rand_distr/tests/cdf.rs @@ -434,11 +434,14 @@ fn poisson() { 1.0 - lambda.inc_gamma(k as f64 + 1.0) } - let parameters = [0.5, 1.0, 7.5, 32.0, 100.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)); + test_discrete(seed as u64, dist, |k| cdf(k, lambda as f64)); } } From 48b6156fa86505844bc35e3d65ca941e1b55d54f Mon Sep 17 00:00:00 2001 From: JamboChen Date: Thu, 10 Oct 2024 17:09:40 +0800 Subject: [PATCH 08/24] add: `triangular` test --- rand_distr/tests/cdf.rs | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/rand_distr/tests/cdf.rs b/rand_distr/tests/cdf.rs index b0ecb9be4b..61c4e92ba2 100644 --- a/rand_distr/tests/cdf.rs +++ b/rand_distr/tests/cdf.rs @@ -338,6 +338,34 @@ fn beta() { } } +#[test] +fn triangular() { + fn cdf(x: f64, a: f64, b: f64, c: f64) -> f64 { + if x <= a { + 0.0 + } else if a < x && x <= c { + (x - a).powi(2) / ((b - a) * (c - a)) + } else if c < x && x < b { + 1.0 - (b - x).powi(2) / ((b - a) * (b - c)) + } else { + 1.0 + } + } + + let parameters = [ + (0.0, 1.0, 0.0001), + (0.0, 1.0, 0.9999), + (0.0, 1.0, 0.5), + (0.0, 100.0, 50.0), + (-100.0, 100.0, 0.0), + ]; + + for (seed, (a, b, c)) in parameters.into_iter().enumerate() { + let dist = rand_distr::Triangular::new(a, b, c).unwrap(); + test_continuous(seed as u64, dist, |x| cdf(x, a, b, c)); + } +} + fn binomial_cdf(k: i64, p: f64, n: u64) -> f64 { if k < 0 { return 0.0; From 10c19b0059b81c6a4b8f004e3dc581beb8532f71 Mon Sep 17 00:00:00 2001 From: JamboChen Date: Thu, 10 Oct 2024 21:15:36 +0800 Subject: [PATCH 09/24] Refactor: poisson function to use `gamma_lr` --- rand_distr/tests/cdf.rs | 79 +++++++++++++++++++++++++++++++++++++++-- 1 file changed, 76 insertions(+), 3 deletions(-) diff --git a/rand_distr/tests/cdf.rs b/rand_distr/tests/cdf.rs index 61c4e92ba2..3f95a1b100 100644 --- a/rand_distr/tests/cdf.rs +++ b/rand_distr/tests/cdf.rs @@ -459,11 +459,11 @@ fn poisson() { return 0.0; } - 1.0 - lambda.inc_gamma(k as f64 + 1.0) + 1.0 - gamma_lr(k as f64 + 1.0, lambda) } - let parameters = [ - 0.1_f32, 1.0, 7.5, + 0.1_f32, 1.0, 7.5, 2e4, + // 1e9, // 1.844E+19, // fail case ]; @@ -480,3 +480,76 @@ fn ln_factorial(n: u64) -> f64 { fn ln_binomial(n: u64, k: u64) -> f64 { ln_factorial(n) - ln_factorial(k) - ln_factorial(n - k) } + +/// https://docs.rs/statrs/latest/src/statrs/function/gamma.rs.html#260-347 +fn gamma_lr(a: f64, x: f64) -> f64 { + let eps = 0.000000000000001; + let big = 4503599627370496.0; + let big_inv = 2.22044604925031308085e-16; + + let ax = a * x.ln() - x - a.lgamma().0; + if ax < -709.78271289338399 { + if a < x { + return 1.0; + } + return 0.0; + } + if x <= 1.0 || x <= a { + let mut r2 = a; + let mut c2 = 1.0; + let mut ans2 = 1.0; + loop { + r2 += 1.0; + c2 *= x / r2; + ans2 += c2; + + if c2 / ans2 <= eps { + break; + } + } + return ax.exp() * ans2 / a; + } + + let mut y = 1.0 - a; + let mut z = x + y + 1.0; + let mut c = 0; + + let mut p3 = 1.0; + let mut q3 = x; + let mut p2 = x + 1.0; + let mut q2 = z * x; + let mut ans = p2 / q2; + + loop { + y += 1.0; + z += 2.0; + c += 1; + let yc = y * f64::from(c); + + let p = p2 * z - p3 * yc; + let q = q2 * z - q3 * yc; + + p3 = p2; + p2 = p; + q3 = q2; + q2 = q; + + if p.abs() > big { + p3 *= big_inv; + p2 *= big_inv; + q3 *= big_inv; + q2 *= big_inv; + } + + if q != 0.0 { + let nextans = p / q; + let error = ((ans - nextans) / nextans).abs(); + ans = nextans; + + if error <= eps { + break; + } + } + } + 1.0 - ax.exp() * ans +} From d22304d7ae015aafeb43df8b914e90de56742c3d Mon Sep 17 00:00:00 2001 From: JamboChen Date: Thu, 10 Oct 2024 21:34:34 +0800 Subject: [PATCH 10/24] add: ks test for `Pareto` --- rand_distr/tests/cdf.rs | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/rand_distr/tests/cdf.rs b/rand_distr/tests/cdf.rs index 3f95a1b100..2edf42bd48 100644 --- a/rand_distr/tests/cdf.rs +++ b/rand_distr/tests/cdf.rs @@ -227,6 +227,32 @@ fn log_normal() { } } +#[test] +fn pareto() { + fn cdf(x: f64, scale: f64, alpha: f64) -> f64 { + if x <= scale { + 0.0 + } else { + 1.0 - (scale / x).powf(alpha) + } + } + + let parameters = [ + (1.0, 1.0), + (1.0, 0.1), + (1.0, 10.0), + (1.0, 100.0), + (0.1, 1.0), + (10.0, 1.0), + (100.0, 1.0), + ]; + + for (seed, (scale, alpha)) in parameters.into_iter().enumerate() { + let dist = rand_distr::Pareto::new(scale, alpha).unwrap(); + test_continuous(seed as u64, dist, |x| cdf(x, scale, alpha)); + } +} + #[test] fn exp() { fn cdf(x: f64, lambda: f64) -> f64 { From c71b3ffe876e10b92fc614d8e4279aa7da839b89 Mon Sep 17 00:00:00 2001 From: JamboChen Date: Thu, 10 Oct 2024 21:46:59 +0800 Subject: [PATCH 11/24] add: ks test for `Gumbel` --- rand_distr/tests/cdf.rs | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/rand_distr/tests/cdf.rs b/rand_distr/tests/cdf.rs index 2edf42bd48..031037542c 100644 --- a/rand_distr/tests/cdf.rs +++ b/rand_distr/tests/cdf.rs @@ -292,6 +292,26 @@ fn weibull() { } } +#[test] +fn gumbel() { + fn cdf(x: f64, mu: f64, beta: f64) -> f64 { + (-(-(x - mu) / beta).exp()).exp() + } + + let parameters = [ + (0.0, 1.0), + (1.0, 2.0), + (-1.0, 0.5), + (10.0, 0.1), + (100.0, 0.0001), + ]; + + for (seed, (mu, beta)) in parameters.into_iter().enumerate() { + let dist = rand_distr::Gumbel::new(mu, beta).unwrap(); + test_continuous(seed as u64, dist, |x| cdf(x, mu, beta)); + } +} + #[test] fn gamma() { fn cdf(x: f64, shape: f64, scale: f64) -> f64 { From 02d05e92d393c431a576dc2f600a02fc9569f63c Mon Sep 17 00:00:00 2001 From: JamboChen Date: Fri, 11 Oct 2024 00:39:08 +0800 Subject: [PATCH 12/24] add: ks test for `Frechet`, `t`, `F` --- rand_distr/tests/cdf.rs | 66 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 66 insertions(+) diff --git a/rand_distr/tests/cdf.rs b/rand_distr/tests/cdf.rs index 031037542c..df568a5d89 100644 --- a/rand_distr/tests/cdf.rs +++ b/rand_distr/tests/cdf.rs @@ -312,6 +312,30 @@ fn gumbel() { } } +#[test] +fn frechet() { + fn cdf(x: f64, alpha: f64, s: f64, m: f64) -> f64 { + if x < m { + return 0.0; + } + + (-((x - m) / s).powf(-alpha)).exp() + } + + let parameters = [ + (0.5, 2.0, 1.0), + (1.0, 1.0, 1.0), + (10.0, 0.1, 1.0), + (100.0, 0.0001, 1.0), + (0.9999, 2.0, 1.0), + ]; + + for (seed, (alpha, s, m)) in parameters.into_iter().enumerate() { + let dist = rand_distr::Frechet::new(m, s, alpha).unwrap(); + test_continuous(seed as u64, dist, |x| cdf(x, alpha, s, m)); + } +} + #[test] fn gamma() { fn cdf(x: f64, shape: f64, scale: f64) -> f64 { @@ -356,6 +380,48 @@ fn chi_squared() { test_continuous(seed as u64, dist, |x| cdf(x, k)); } } +#[test] +fn studend_t() { + fn cdf(x: f64, df: f64) -> f64 { + let h = df / (df + x.powi(2)); + let ib = 0.5 * h.inc_beta(df / 2.0, 0.5, 0.5.ln_beta(df / 2.0)); + if x < 0.0 { + ib + } else { + 1.0 - ib + } + } + + let parameters = [1.0, 10.0, 50.0]; + + for (seed, df) in parameters.into_iter().enumerate() { + let dist = rand_distr::StudentT::new(df).unwrap(); + test_continuous(seed as u64, dist, |x| cdf(x, df)); + } +} + +#[test] +fn fisher_f() { + fn cdf(x: f64, m: f64, n: f64) -> f64 { + if m == 1.0 && x <= 0.0 { + return 0.0; + } else if x < 0.0 { + return 0.0; + } else { + let k = m * x / (m * x + n); + let d1 = m / 2.0; + let d2 = n / 2.0; + k.inc_beta(d1, d2, d1.ln_beta(d2)) + } + } + + let parameters = [(1.0, 1.0), (1.0, 2.0), (2.0, 1.0), (50.0, 1.0)]; + + for (seed, (m, n)) in parameters.into_iter().enumerate() { + let dist = rand_distr::FisherF::new(m, n).unwrap(); + test_continuous(seed as u64, dist, |x| cdf(x, m, n)); + } +} #[test] fn beta() { From 2bf4dad0ec4a0c2ef9bf023252ef34886b118bd9 Mon Sep 17 00:00:00 2001 From: JamboChen Date: Fri, 11 Oct 2024 00:42:05 +0800 Subject: [PATCH 13/24] clippy --- rand_distr/tests/cdf.rs | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/rand_distr/tests/cdf.rs b/rand_distr/tests/cdf.rs index df568a5d89..fb8f3ffd63 100644 --- a/rand_distr/tests/cdf.rs +++ b/rand_distr/tests/cdf.rs @@ -403,10 +403,8 @@ fn studend_t() { #[test] fn fisher_f() { fn cdf(x: f64, m: f64, n: f64) -> f64 { - if m == 1.0 && x <= 0.0 { - return 0.0; - } else if x < 0.0 { - return 0.0; + if (m == 1.0 && x <= 0.0) || x < 0.0 { + 0.0 } else { let k = m * x / (m * x + n); let d1 = m / 2.0; @@ -597,10 +595,10 @@ fn ln_binomial(n: u64, k: u64) -> f64 { fn gamma_lr(a: f64, x: f64) -> f64 { let eps = 0.000000000000001; let big = 4503599627370496.0; - let big_inv = 2.22044604925031308085e-16; + let big_inv = 2.220_446_049_250_313e-16; let ax = a * x.ln() - x - a.lgamma().0; - if ax < -709.78271289338399 { + if ax < -709.782_712_893_384 { if a < x { return 1.0; } From 046765531b62ef2fc078f14c2c6b07278d083e5d Mon Sep 17 00:00:00 2001 From: JamboChen Date: Fri, 11 Oct 2024 16:45:26 +0800 Subject: [PATCH 14/24] add: ks test for `zeta`, `zipf` --- rand_distr/tests/cdf.rs | 85 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 85 insertions(+) diff --git a/rand_distr/tests/cdf.rs b/rand_distr/tests/cdf.rs index fb8f3ffd63..7951323273 100644 --- a/rand_distr/tests/cdf.rs +++ b/rand_distr/tests/cdf.rs @@ -336,6 +336,46 @@ fn frechet() { } } +#[test] +fn zeta() { + fn cdf(k: i64, s: f64) -> f64 { + if k <= 1 { + return 0.0; + } + + gen_harmonic(k as u64, s) / zeta_func(s) + } + + let parameters = [ + // 2.0, 3.7, 5.0, 100.0, // all fail + ]; + + for (seed, s) in parameters.into_iter().enumerate() { + let dist = rand_distr::Zeta::new(s).unwrap(); + test_discrete(seed as u64, dist, |k| cdf(k, s)); + } +} + +#[test] +fn zipf() { + fn cdf(k: i64, n: u64, s: f64) -> f64 { + if k < 1 { + return 0.0; + } + if k > n as i64 { + return 1.0; + } + gen_harmonic(k as u64, s) / gen_harmonic(n, s) + } + + let parameters = [(1000, 1.0), (500, 2.0), (1000, 0.5)]; + + for (seed, (n, x)) in parameters.into_iter().enumerate() { + let dist = rand_distr::Zipf::new(n, x).unwrap(); + test_discrete(seed as u64, dist, |k| cdf(k, n, x)); + } +} + #[test] fn gamma() { fn cdf(x: f64, shape: f64, scale: f64) -> f64 { @@ -663,3 +703,48 @@ fn gamma_lr(a: f64, x: f64) -> f64 { } 1.0 - ax.exp() * ans } + +fn gen_harmonic(n: u64, m: f64) -> f64 { + match n { + 0 => 1.0, + _ => (0..n).fold(0.0, |acc, x| acc + (x as f64 + 1.0).powf(-m)), + } +} + +/// https://docs.rs/spfunc/latest/src/spfunc/zeta.rs.html#20-43 +fn zeta_func(s: f64) -> f64 { + if s == 1.0 { + return f64::INFINITY; + } + + let mut akn = vec![1.0]; + let mut two_pow = 0.5; + let head = 1.0 / (1.0 - 2.0.powf(1.0 - s)); + let mut tail_prev = 0.0; + let mut tail = two_pow * akn[0]; + + while (tail - tail_prev).abs() >= f64::EPSILON { + update_akn(&mut akn, s); + two_pow /= 2.0; + tail_prev = tail; + tail += two_pow * akn.iter().sum::(); + } + + head * tail +} + +fn update_akn(akn: &mut Vec, s: f64) { + let n = akn.len() - 1; + + let n1 = n as f64 + 1.0; + + akn.iter_mut().enumerate().for_each(|(k, a)| { + let num = n1; + let den = n + 1 - k; + *a *= num / den as f64; + }); + let p1 = -1.0 / n1; + let p2 = (n1 / (n1 + 1.0)).powf(s); + + akn.push(p1 * p2 * akn[n]); +} From e67f8022db615572b35a16c00d17cf35e1c2621b Mon Sep 17 00:00:00 2001 From: JamboChen Date: Fri, 11 Oct 2024 17:09:57 +0800 Subject: [PATCH 15/24] fix: zeta cdf's support --- rand_distr/tests/cdf.rs | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/rand_distr/tests/cdf.rs b/rand_distr/tests/cdf.rs index 7951323273..3203a6a7a6 100644 --- a/rand_distr/tests/cdf.rs +++ b/rand_distr/tests/cdf.rs @@ -339,16 +339,14 @@ fn frechet() { #[test] fn zeta() { fn cdf(k: i64, s: f64) -> f64 { - if k <= 1 { + if k < 1 { return 0.0; } gen_harmonic(k as u64, s) / zeta_func(s) } - let parameters = [ - // 2.0, 3.7, 5.0, 100.0, // all fail - ]; + let parameters = [2.0, 3.7, 5.0, 100.0]; for (seed, s) in parameters.into_iter().enumerate() { let dist = rand_distr::Zeta::new(s).unwrap(); From a8e3596812cfc59ecdb73cebdda7294c59854e12 Mon Sep 17 00:00:00 2001 From: JamboChen Date: Fri, 11 Oct 2024 22:38:53 +0800 Subject: [PATCH 16/24] add: ks test for `SkewNormal` --- rand_distr/tests/cdf.rs | 256 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 256 insertions(+) diff --git a/rand_distr/tests/cdf.rs b/rand_distr/tests/cdf.rs index 3203a6a7a6..a7801684df 100644 --- a/rand_distr/tests/cdf.rs +++ b/rand_distr/tests/cdf.rs @@ -158,6 +158,21 @@ fn normal() { } } +#[test] +fn skew_normal() { + fn cdf(x: f64, location: f64, scale: f64, shape: f64) -> f64 { + let norm = (x - location) / scale; + phi(norm) - 2.0 * owen_t(norm, shape) + } + + let parameters = [(0.0, 1.0, 5.0), (1.0, 10.0, -5.0), (-1.0, 0.00001, 0.0)]; + + for (seed, (location, scale, shape)) in parameters.into_iter().enumerate() { + let dist = rand_distr::SkewNormal::new(location, scale, shape).unwrap(); + test_continuous(seed as u64, dist, |x| cdf(x, location, scale, shape)); + } +} + #[test] fn cauchy() { fn cdf(x: f64, median: f64, scale: f64) -> f64 { @@ -746,3 +761,244 @@ fn update_akn(akn: &mut Vec, s: f64) { akn.push(p1 * p2 * akn[n]); } + + +/// [1] FAST AND ACCURATE CALCULATION OF OWEN’S T-FUNCTION +/// by Mike Patefield and David Tandy +/// http://www.jstatsoft.org/v05/i05/paper + +fn owen_t(h: f64, a: f64) -> f64 { + let absh = h.abs(); + let absa = a.abs(); + let ah = absa * absh; + + let mut t; + if absa <= 1.0 { + t = tf(absh, absa, ah); + } else if absh <= 0.67 { + t = 0.25 - znorm1(absh) * znorm1(ah) - tf(ah, 1.0 / absa, absh); + } else { + let normh = znorm2(absh); + let normah = znorm2(ah); + t = 0.5 * (normh + normah) - normh * normah - tf(ah, 1.0 / absa, absh); + } + + if a < 0.0 { + t = -t; + } + + t +} + +fn tf(h: f64, a: f64, ah: f64) -> f64 { + let rtwopi = 0.15915494309189533577; + let rrtpi = 0.39894228040143267794; + + let c2 = [ + 0.99999999999999987510, + -0.99999999999988796462, + 0.99999999998290743652, + -0.99999999896282500134, + 0.99999996660459362918, + -0.99999933986272476760, + 0.99999125611136965852, + -0.99991777624463387686, + 0.99942835555870132569, + -0.99697311720723000295, + 0.98751448037275303682, + -0.95915857980572882813, + 0.89246305511006708555, + -0.76893425990463999675, + 0.58893528468484693250, + -0.38380345160440256652, + 0.20317601701045299653, + -0.82813631607004984866e-01, + 0.24167984735759576523e-01, + -0.44676566663971825242e-02, + 0.39141169402373836468e-03, + ]; + + let pts = [ + 0.35082039676451715489e-02, + 0.31279042338030753740e-01, + 0.85266826283219451090e-01, + 0.16245071730812277011, + 0.25851196049125434828, + 0.36807553840697533536, + 0.48501092905604697475, + 0.60277514152618576821, + 0.71477884217753226516, + 0.81475510988760098605, + 0.89711029755948965867, + 0.95723808085944261843, + 0.99178832974629703586, + ]; + + let wts = [ + 0.18831438115323502887e-01, + 0.18567086243977649478e-01, + 0.18042093461223385584e-01, + 0.17263829606398753364e-01, + 0.16243219975989856730e-01, + 0.14994592034116704829e-01, + 0.13535474469662088392e-01, + 0.11886351605820165233e-01, + 0.10070377242777431897e-01, + 0.81130545742299586629e-02, + 0.60419009528470238773e-02, + 0.38862217010742057883e-02, + 0.16793031084546090448e-02, + ]; + + let hrange = [ + 0.02, 0.06, 0.09, 0.125, 0.26, 0.4, 0.6, 1.6, 1.7, 2.33, 2.4, 3.36, 3.4, 4.8, + ]; + let arange = [0.025, 0.09, 0.15, 0.36, 0.5, 0.9, 0.99999]; + + let select = [ + [1, 1, 2, 13, 13, 13, 13, 13, 13, 13, 13, 16, 16, 16, 9], + [1, 2, 2, 3, 3, 5, 5, 14, 14, 15, 15, 16, 16, 16, 9], + [2, 2, 3, 3, 3, 5, 5, 15, 15, 15, 15, 16, 16, 16, 10], + [2, 2, 3, 5, 5, 5, 5, 7, 7, 16, 16, 16, 16, 16, 10], + [2, 3, 3, 5, 5, 6, 6, 8, 8, 17, 17, 17, 12, 12, 11], + [2, 3, 5, 5, 5, 6, 6, 8, 8, 17, 17, 17, 12, 12, 12], + [2, 3, 4, 4, 6, 6, 8, 8, 17, 17, 17, 17, 17, 12, 12], + [2, 3, 4, 4, 6, 6, 18, 18, 18, 18, 17, 17, 17, 12, 12], + ]; + + let ihint = match hrange.iter().position(|&r| h < r) { + Some(idx) => idx, + None => 14, + }; + + let iaint = match arange.iter().position(|&r| a < r) { + Some(idx) => idx, + None => 7, + }; + + let icode = select[iaint][ihint]; + let m = [ + 2, 3, 4, 5, 7, 10, 12, 18, 10, 20, 30, 20, 4, 7, 8, 20, 13, 0, + ][icode - 1]; + let method = [1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 4, 4, 4, 4, 5, 6][icode - 1]; + + match method { + 1 => { + let hs = -0.5 * h * h; + let dhs = hs.exp(); + let as_ = a * a; + let mut j = 1; + let mut jj = 1; + let mut aj = rtwopi * a; + let mut tf = rtwopi * a.atan(); + let mut dj = dhs - 1.0; + let mut gj = hs * dhs; + loop { + tf += dj * aj / (jj as f64); + if j >= m { + return tf; + } + j += 1; + jj += 2; + aj *= as_; + dj = gj - dj; + gj *= hs / (j as f64); + } + } + 2 => { + let maxii = m + m + 1; + let mut ii = 1; + let mut tf = 0.0; + let hs = h * h; + let as_ = -a * a; + let mut vi = rrtpi * a * (-0.5 * ah * ah).exp(); + let mut z = znorm1(ah) / h; + let y = 1.0 / hs; + loop { + tf += z; + if ii >= maxii { + tf *= rrtpi * (-0.5 * hs).exp(); + return tf; + } + z = y * (vi - (ii as f64) * z); + vi *= as_; + ii += 2; + } + } + 3 => { + let mut i = 1; + let mut ii = 1; + let mut tf = 0.0; + let hs = h * h; + let as_ = a * a; + let mut vi = rrtpi * a * (-0.5 * ah * ah).exp(); + let mut zi = znorm1(ah) / h; + let y = 1.0 / hs; + loop { + tf += zi * c2[i - 1]; + if i > m { + tf *= rrtpi * (-0.5 * hs).exp(); + return tf; + } + zi = y * ((ii as f64) * zi - vi); + vi *= as_; + i += 1; + ii += 2; + } + } + 4 => { + let maxii = m + m + 1; + let mut ii = 1; + let mut tf = 0.0; + let hs = h * h; + let as_ = -a * a; + let mut ai = rtwopi * a * (-0.5 * hs * (1.0 - as_)).exp(); + let mut yi = 1.0; + loop { + tf += ai * yi; + if ii >= maxii { + return tf; + } + ii += 2; + yi = (1.0 - hs * yi) / (ii as f64); + ai *= as_; + } + } + 5 => { + let mut tf = 0.0; + let as_ = a * a; + let hs = -0.5 * h * h; + for i in 0..m { + let r = 1.0 + as_ * pts[i]; + tf += wts[i] * (hs * r).exp() / r; + } + tf *= a; + return tf; + } + 6 => { + let normh = znorm2(h); + let mut tf = 0.5 * normh * (1.0 - normh); + let y = 1.0 - a; + let r = (y / (1.0 + a)).atan(); + if r != 0.0 { + tf -= rtwopi * r * (-0.5 * y * h * h / r).exp(); + } + return tf; + } + _ => 0.0, + } +} + +fn phi(x: f64) -> f64 { + normal_cdf(x, 0.0, 1.0) +} + +// P(0 ≤ Z ≤ x) +fn znorm1(x: f64) -> f64 { + phi(x) - 0.5 +} + +// P(x ≤ Z < ∞) +fn znorm2(x: f64) -> f64 { + 1.0 - phi(x) +} From 88e58bca8c4f01f8b0e9371b382a2e0c7f0aade8 Mon Sep 17 00:00:00 2001 From: JamboChen Date: Fri, 11 Oct 2024 22:39:48 +0800 Subject: [PATCH 17/24] clippy --- rand_distr/tests/cdf.rs | 112 +++++++++++++++++++--------------------- 1 file changed, 53 insertions(+), 59 deletions(-) diff --git a/rand_distr/tests/cdf.rs b/rand_distr/tests/cdf.rs index a7801684df..aeedfec0ba 100644 --- a/rand_distr/tests/cdf.rs +++ b/rand_distr/tests/cdf.rs @@ -791,63 +791,63 @@ fn owen_t(h: f64, a: f64) -> f64 { } fn tf(h: f64, a: f64, ah: f64) -> f64 { - let rtwopi = 0.15915494309189533577; - let rrtpi = 0.39894228040143267794; + let rtwopi = 0.159_154_943_091_895_35; + let rrtpi = 0.398_942_280_401_432_7; let c2 = [ - 0.99999999999999987510, - -0.99999999999988796462, - 0.99999999998290743652, - -0.99999999896282500134, - 0.99999996660459362918, - -0.99999933986272476760, - 0.99999125611136965852, - -0.99991777624463387686, - 0.99942835555870132569, - -0.99697311720723000295, - 0.98751448037275303682, - -0.95915857980572882813, - 0.89246305511006708555, - -0.76893425990463999675, - 0.58893528468484693250, - -0.38380345160440256652, - 0.20317601701045299653, - -0.82813631607004984866e-01, - 0.24167984735759576523e-01, - -0.44676566663971825242e-02, - 0.39141169402373836468e-03, + 0.999_999_999_999_999_9, + -0.999_999_999_999_888, + 0.999_999_999_982_907_5, + -0.999_999_998_962_825, + 0.999_999_966_604_593_7, + -0.999_999_339_862_724_7, + 0.999_991_256_111_369_6, + -0.999_917_776_244_633_8, + 0.999_428_355_558_701_4, + -0.996_973_117_207_23, + 0.987_514_480_372_753, + -0.959_158_579_805_728_8, + 0.892_463_055_110_067_1, + -0.768_934_259_904_64, + 0.588_935_284_684_846_9, + -0.383_803_451_604_402_55, + 0.203_176_017_010_453, + -8.281_363_160_700_499e-2, + 2.416_798_473_575_957_8e-2, + -4.467_656_666_397_183e-3, + 3.914_116_940_237_383_6e-4, ]; let pts = [ - 0.35082039676451715489e-02, - 0.31279042338030753740e-01, - 0.85266826283219451090e-01, - 0.16245071730812277011, - 0.25851196049125434828, - 0.36807553840697533536, - 0.48501092905604697475, - 0.60277514152618576821, - 0.71477884217753226516, - 0.81475510988760098605, - 0.89711029755948965867, - 0.95723808085944261843, - 0.99178832974629703586, + 3.508_203_967_645_171_6e-3, + 3.127_904_233_803_075_6e-2, + 8.526_682_628_321_945e-2, + 0.162_450_717_308_122_77, + 0.258_511_960_491_254_36, + 0.368_075_538_406_975_3, + 0.485_010_929_056_047, + 0.602_775_141_526_185_7, + 0.714_778_842_177_532_3, + 0.814_755_109_887_601, + 0.897_110_297_559_489_7, + 0.957_238_080_859_442_6, + 0.991_788_329_746_297, ]; let wts = [ - 0.18831438115323502887e-01, - 0.18567086243977649478e-01, - 0.18042093461223385584e-01, - 0.17263829606398753364e-01, - 0.16243219975989856730e-01, - 0.14994592034116704829e-01, - 0.13535474469662088392e-01, - 0.11886351605820165233e-01, - 0.10070377242777431897e-01, - 0.81130545742299586629e-02, - 0.60419009528470238773e-02, - 0.38862217010742057883e-02, - 0.16793031084546090448e-02, + 1.883_143_811_532_350_3e-2, + 1.856_708_624_397_765e-2, + 1.804_209_346_122_338_5e-2, + 1.726_382_960_639_875_2e-2, + 1.624_321_997_598_985_8e-2, + 1.499_459_203_411_670_5e-2, + 1.353_547_446_966_209e-2, + 1.188_635_160_582_016_5e-2, + 1.007_037_724_277_743_2e-2, + 8.113_054_574_229_958e-3, + 6.041_900_952_847_024e-3, + 3.886_221_701_074_205_7e-3, + 1.679_303_108_454_609e-3, ]; let hrange = [ @@ -866,15 +866,9 @@ fn tf(h: f64, a: f64, ah: f64) -> f64 { [2, 3, 4, 4, 6, 6, 18, 18, 18, 18, 17, 17, 17, 12, 12], ]; - let ihint = match hrange.iter().position(|&r| h < r) { - Some(idx) => idx, - None => 14, - }; + let ihint = hrange.iter().position(|&r| h < r).unwrap_or(14); - let iaint = match arange.iter().position(|&r| a < r) { - Some(idx) => idx, - None => 7, - }; + let iaint = arange.iter().position(|&r| a < r).unwrap_or(7); let icode = select[iaint][ihint]; let m = [ @@ -973,7 +967,7 @@ fn tf(h: f64, a: f64, ah: f64) -> f64 { tf += wts[i] * (hs * r).exp() / r; } tf *= a; - return tf; + tf } 6 => { let normh = znorm2(h); @@ -983,7 +977,7 @@ fn tf(h: f64, a: f64, ah: f64) -> f64 { if r != 0.0 { tf -= rtwopi * r * (-0.5 * y * h * h / r).exp(); } - return tf; + tf } _ => 0.0, } From dbfda2ba5644cf06b7220f9c9be0b3db1acd9ba8 Mon Sep 17 00:00:00 2001 From: JamboChen Date: Fri, 11 Oct 2024 22:41:21 +0800 Subject: [PATCH 18/24] fmt --- rand_distr/tests/cdf.rs | 1 - 1 file changed, 1 deletion(-) diff --git a/rand_distr/tests/cdf.rs b/rand_distr/tests/cdf.rs index aeedfec0ba..053f63daad 100644 --- a/rand_distr/tests/cdf.rs +++ b/rand_distr/tests/cdf.rs @@ -762,7 +762,6 @@ fn update_akn(akn: &mut Vec, s: f64) { akn.push(p1 * p2 * akn[n]); } - /// [1] FAST AND ACCURATE CALCULATION OF OWEN’S T-FUNCTION /// by Mike Patefield and David Tandy /// http://www.jstatsoft.org/v05/i05/paper From 102411d5d0f3ff43a444796b45e14182719c585b Mon Sep 17 00:00:00 2001 From: JamboChen Date: Thu, 17 Oct 2024 19:40:50 +0800 Subject: [PATCH 19/24] Delete unreasonable test parameters according to comments --- rand_distr/tests/cdf.rs | 17 +++-------------- 1 file changed, 3 insertions(+), 14 deletions(-) diff --git a/rand_distr/tests/cdf.rs b/rand_distr/tests/cdf.rs index 053f63daad..d334dd8021 100644 --- a/rand_distr/tests/cdf.rs +++ b/rand_distr/tests/cdf.rs @@ -234,8 +234,6 @@ fn log_normal() { (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)); @@ -298,7 +296,7 @@ fn weibull() { (10.0, 0.1), (0.1, 10.0), (15.0, 20.0), - // (1000.0, 0.001), // Fail case + (1000.0, 0.01), ]; for (seed, (lambda, k)) in parameters.into_iter().enumerate() { @@ -423,10 +421,7 @@ fn chi_squared() { (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, - ]; + let parameters = [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(); @@ -487,13 +482,7 @@ fn 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 - ]; + let parameters = [(0.5, 0.5), (2.0, 3.5), (10.0, 1.0), (100.0, 50.0)]; for (seed, (alpha, beta)) in parameters.into_iter().enumerate() { let dist = rand_distr::Beta::new(alpha, beta).unwrap(); From ac7e19a5ed71a4a6fa57fb7f59af3830123918ae Mon Sep 17 00:00:00 2001 From: JamboChen Date: Thu, 17 Oct 2024 21:22:45 +0800 Subject: [PATCH 20/24] fix: `f64` poisson distribution on 1.16.1 --- rand_distr/tests/cdf.rs | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/rand_distr/tests/cdf.rs b/rand_distr/tests/cdf.rs index d334dd8021..b387d5bc98 100644 --- a/rand_distr/tests/cdf.rs +++ b/rand_distr/tests/cdf.rs @@ -120,11 +120,12 @@ pub fn test_continuous(seed: u64, dist: impl Distribution, cdf: impl Fn(f64 /// Tests a distribution over integers against an analytical CDF. /// The analytical CDF must not have jump points which are not integers. -pub fn test_discrete>( - seed: u64, - dist: impl Distribution, - cdf: impl Fn(i64) -> f64, -) { +pub fn test_discrete(seed: u64, dist: D, cdf: F) +where + I: AsPrimitive, + D: Distribution, + F: Fn(i64) -> f64, +{ let ecdf = sample_ecdf(seed, dist); let ks_statistic = kolmogorov_smirnov_statistic_discrete(ecdf, cdf); @@ -606,6 +607,7 @@ fn hypergeometric() { #[test] fn poisson() { + use rand_distr::Poisson; fn cdf(k: i64, lambda: f64) -> f64 { if k < 0 || lambda <= 0.0 { return 0.0; @@ -614,14 +616,14 @@ fn poisson() { 1.0 - gamma_lr(k as f64 + 1.0, lambda) } let parameters = [ - 0.1_f32, 1.0, 7.5, 2e4, - // 1e9, + 0.1, 1.0, 7.5, + // 1e9, passed case but too slow // 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)); + let dist = Poisson::new(lambda).unwrap(); + test_discrete::, _>(seed as u64, dist, |k| cdf(k, lambda)); } } From 436ab975f91fcceb652b6eb4d1fe02a5651aa9f7 Mon Sep 17 00:00:00 2001 From: JamboChen Date: Thu, 17 Oct 2024 23:11:52 +0800 Subject: [PATCH 21/24] Add: `CHANGELOG.md` entry --- rand_distr/CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/rand_distr/CHANGELOG.md b/rand_distr/CHANGELOG.md index a19641752f..c9ab729be8 100644 --- a/rand_distr/CHANGELOG.md +++ b/rand_distr/CHANGELOG.md @@ -9,6 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Fix panic in Binomial (#1484) - Move some of the computations in Binomial from `sample` to `new` (#1484) - Add Kolmogorov Smirnov test for sampling of `Normal` and `Binomial` (#1494) +- Add Kolmogorov Smirnov test for more distributions (#1504) ### Added - Add plots for `rand_distr` distributions to documentation (#1434) From 771e79ba3255e791b5463dcb518cb9b39528a0ee Mon Sep 17 00:00:00 2001 From: JamboChen Date: Fri, 18 Oct 2024 00:36:53 +0800 Subject: [PATCH 22/24] Tidy up `owen_t` functions --- rand_distr/tests/cdf.rs | 394 ++++++++++++++++++++-------------------- 1 file changed, 198 insertions(+), 196 deletions(-) diff --git a/rand_distr/tests/cdf.rs b/rand_distr/tests/cdf.rs index b387d5bc98..fde345a124 100644 --- a/rand_distr/tests/cdf.rs +++ b/rand_distr/tests/cdf.rs @@ -753,10 +753,11 @@ fn update_akn(akn: &mut Vec, s: f64) { akn.push(p1 * p2 * akn[n]); } -/// [1] FAST AND ACCURATE CALCULATION OF OWEN’S T-FUNCTION -/// by Mike Patefield and David Tandy -/// http://www.jstatsoft.org/v05/i05/paper - +/// [1] Patefield, M. (2000). Fast and Accurate Calculation of Owen’s T Function. +/// Journal of Statistical Software, 5(5), 1–25. +/// https://doi.org/10.18637/jss.v005.i05 +/// +/// This function is ported to Rust from the Fortran code provided in the paper fn owen_t(h: f64, a: f64) -> f64 { let absh = h.abs(); let absa = a.abs(); @@ -777,212 +778,213 @@ fn owen_t(h: f64, a: f64) -> f64 { t = -t; } - t -} - -fn tf(h: f64, a: f64, ah: f64) -> f64 { - let rtwopi = 0.159_154_943_091_895_35; - let rrtpi = 0.398_942_280_401_432_7; - - let c2 = [ - 0.999_999_999_999_999_9, - -0.999_999_999_999_888, - 0.999_999_999_982_907_5, - -0.999_999_998_962_825, - 0.999_999_966_604_593_7, - -0.999_999_339_862_724_7, - 0.999_991_256_111_369_6, - -0.999_917_776_244_633_8, - 0.999_428_355_558_701_4, - -0.996_973_117_207_23, - 0.987_514_480_372_753, - -0.959_158_579_805_728_8, - 0.892_463_055_110_067_1, - -0.768_934_259_904_64, - 0.588_935_284_684_846_9, - -0.383_803_451_604_402_55, - 0.203_176_017_010_453, - -8.281_363_160_700_499e-2, - 2.416_798_473_575_957_8e-2, - -4.467_656_666_397_183e-3, - 3.914_116_940_237_383_6e-4, - ]; - - let pts = [ - 3.508_203_967_645_171_6e-3, - 3.127_904_233_803_075_6e-2, - 8.526_682_628_321_945e-2, - 0.162_450_717_308_122_77, - 0.258_511_960_491_254_36, - 0.368_075_538_406_975_3, - 0.485_010_929_056_047, - 0.602_775_141_526_185_7, - 0.714_778_842_177_532_3, - 0.814_755_109_887_601, - 0.897_110_297_559_489_7, - 0.957_238_080_859_442_6, - 0.991_788_329_746_297, - ]; - - let wts = [ - 1.883_143_811_532_350_3e-2, - 1.856_708_624_397_765e-2, - 1.804_209_346_122_338_5e-2, - 1.726_382_960_639_875_2e-2, - 1.624_321_997_598_985_8e-2, - 1.499_459_203_411_670_5e-2, - 1.353_547_446_966_209e-2, - 1.188_635_160_582_016_5e-2, - 1.007_037_724_277_743_2e-2, - 8.113_054_574_229_958e-3, - 6.041_900_952_847_024e-3, - 3.886_221_701_074_205_7e-3, - 1.679_303_108_454_609e-3, - ]; - - let hrange = [ - 0.02, 0.06, 0.09, 0.125, 0.26, 0.4, 0.6, 1.6, 1.7, 2.33, 2.4, 3.36, 3.4, 4.8, - ]; - let arange = [0.025, 0.09, 0.15, 0.36, 0.5, 0.9, 0.99999]; - - let select = [ - [1, 1, 2, 13, 13, 13, 13, 13, 13, 13, 13, 16, 16, 16, 9], - [1, 2, 2, 3, 3, 5, 5, 14, 14, 15, 15, 16, 16, 16, 9], - [2, 2, 3, 3, 3, 5, 5, 15, 15, 15, 15, 16, 16, 16, 10], - [2, 2, 3, 5, 5, 5, 5, 7, 7, 16, 16, 16, 16, 16, 10], - [2, 3, 3, 5, 5, 6, 6, 8, 8, 17, 17, 17, 12, 12, 11], - [2, 3, 5, 5, 5, 6, 6, 8, 8, 17, 17, 17, 12, 12, 12], - [2, 3, 4, 4, 6, 6, 8, 8, 17, 17, 17, 17, 17, 12, 12], - [2, 3, 4, 4, 6, 6, 18, 18, 18, 18, 17, 17, 17, 12, 12], - ]; - - let ihint = hrange.iter().position(|&r| h < r).unwrap_or(14); - - let iaint = arange.iter().position(|&r| a < r).unwrap_or(7); - - let icode = select[iaint][ihint]; - let m = [ - 2, 3, 4, 5, 7, 10, 12, 18, 10, 20, 30, 20, 4, 7, 8, 20, 13, 0, - ][icode - 1]; - let method = [1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 4, 4, 4, 4, 5, 6][icode - 1]; - - match method { - 1 => { - let hs = -0.5 * h * h; - let dhs = hs.exp(); - let as_ = a * a; - let mut j = 1; - let mut jj = 1; - let mut aj = rtwopi * a; - let mut tf = rtwopi * a.atan(); - let mut dj = dhs - 1.0; - let mut gj = hs * dhs; - loop { - tf += dj * aj / (jj as f64); - if j >= m { - return tf; + fn tf(h: f64, a: f64, ah: f64) -> f64 { + let rtwopi = 0.159_154_943_091_895_35; + let rrtpi = 0.398_942_280_401_432_7; + + let c2 = [ + 0.999_999_999_999_999_9, + -0.999_999_999_999_888, + 0.999_999_999_982_907_5, + -0.999_999_998_962_825, + 0.999_999_966_604_593_7, + -0.999_999_339_862_724_7, + 0.999_991_256_111_369_6, + -0.999_917_776_244_633_8, + 0.999_428_355_558_701_4, + -0.996_973_117_207_23, + 0.987_514_480_372_753, + -0.959_158_579_805_728_8, + 0.892_463_055_110_067_1, + -0.768_934_259_904_64, + 0.588_935_284_684_846_9, + -0.383_803_451_604_402_55, + 0.203_176_017_010_453, + -8.281_363_160_700_499e-2, + 2.416_798_473_575_957_8e-2, + -4.467_656_666_397_183e-3, + 3.914_116_940_237_383_6e-4, + ]; + + let pts = [ + 3.508_203_967_645_171_6e-3, + 3.127_904_233_803_075_6e-2, + 8.526_682_628_321_945e-2, + 0.162_450_717_308_122_77, + 0.258_511_960_491_254_36, + 0.368_075_538_406_975_3, + 0.485_010_929_056_047, + 0.602_775_141_526_185_7, + 0.714_778_842_177_532_3, + 0.814_755_109_887_601, + 0.897_110_297_559_489_7, + 0.957_238_080_859_442_6, + 0.991_788_329_746_297, + ]; + + let wts = [ + 1.883_143_811_532_350_3e-2, + 1.856_708_624_397_765e-2, + 1.804_209_346_122_338_5e-2, + 1.726_382_960_639_875_2e-2, + 1.624_321_997_598_985_8e-2, + 1.499_459_203_411_670_5e-2, + 1.353_547_446_966_209e-2, + 1.188_635_160_582_016_5e-2, + 1.007_037_724_277_743_2e-2, + 8.113_054_574_229_958e-3, + 6.041_900_952_847_024e-3, + 3.886_221_701_074_205_7e-3, + 1.679_303_108_454_609e-3, + ]; + + let hrange = [ + 0.02, 0.06, 0.09, 0.125, 0.26, 0.4, 0.6, 1.6, 1.7, 2.33, 2.4, 3.36, 3.4, 4.8, + ]; + let arange = [0.025, 0.09, 0.15, 0.36, 0.5, 0.9, 0.99999]; + + let select = [ + [1, 1, 2, 13, 13, 13, 13, 13, 13, 13, 13, 16, 16, 16, 9], + [1, 2, 2, 3, 3, 5, 5, 14, 14, 15, 15, 16, 16, 16, 9], + [2, 2, 3, 3, 3, 5, 5, 15, 15, 15, 15, 16, 16, 16, 10], + [2, 2, 3, 5, 5, 5, 5, 7, 7, 16, 16, 16, 16, 16, 10], + [2, 3, 3, 5, 5, 6, 6, 8, 8, 17, 17, 17, 12, 12, 11], + [2, 3, 5, 5, 5, 6, 6, 8, 8, 17, 17, 17, 12, 12, 12], + [2, 3, 4, 4, 6, 6, 8, 8, 17, 17, 17, 17, 17, 12, 12], + [2, 3, 4, 4, 6, 6, 18, 18, 18, 18, 17, 17, 17, 12, 12], + ]; + + let ihint = hrange.iter().position(|&r| h < r).unwrap_or(14); + + let iaint = arange.iter().position(|&r| a < r).unwrap_or(7); + + let icode = select[iaint][ihint]; + let m = [ + 2, 3, 4, 5, 7, 10, 12, 18, 10, 20, 30, 20, 4, 7, 8, 20, 13, 0, + ][icode - 1]; + let method = [1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 4, 4, 4, 4, 5, 6][icode - 1]; + + match method { + 1 => { + let hs = -0.5 * h * h; + let dhs = hs.exp(); + let as_ = a * a; + let mut j = 1; + let mut jj = 1; + let mut aj = rtwopi * a; + let mut tf = rtwopi * a.atan(); + let mut dj = dhs - 1.0; + let mut gj = hs * dhs; + loop { + tf += dj * aj / (jj as f64); + if j >= m { + return tf; + } + j += 1; + jj += 2; + aj *= as_; + dj = gj - dj; + gj *= hs / (j as f64); } - j += 1; - jj += 2; - aj *= as_; - dj = gj - dj; - gj *= hs / (j as f64); } - } - 2 => { - let maxii = m + m + 1; - let mut ii = 1; - let mut tf = 0.0; - let hs = h * h; - let as_ = -a * a; - let mut vi = rrtpi * a * (-0.5 * ah * ah).exp(); - let mut z = znorm1(ah) / h; - let y = 1.0 / hs; - loop { - tf += z; - if ii >= maxii { - tf *= rrtpi * (-0.5 * hs).exp(); - return tf; + 2 => { + let maxii = m + m + 1; + let mut ii = 1; + let mut tf = 0.0; + let hs = h * h; + let as_ = -a * a; + let mut vi = rrtpi * a * (-0.5 * ah * ah).exp(); + let mut z = znorm1(ah) / h; + let y = 1.0 / hs; + loop { + tf += z; + if ii >= maxii { + tf *= rrtpi * (-0.5 * hs).exp(); + return tf; + } + z = y * (vi - (ii as f64) * z); + vi *= as_; + ii += 2; } - z = y * (vi - (ii as f64) * z); - vi *= as_; - ii += 2; } - } - 3 => { - let mut i = 1; - let mut ii = 1; - let mut tf = 0.0; - let hs = h * h; - let as_ = a * a; - let mut vi = rrtpi * a * (-0.5 * ah * ah).exp(); - let mut zi = znorm1(ah) / h; - let y = 1.0 / hs; - loop { - tf += zi * c2[i - 1]; - if i > m { - tf *= rrtpi * (-0.5 * hs).exp(); - return tf; + 3 => { + let mut i = 1; + let mut ii = 1; + let mut tf = 0.0; + let hs = h * h; + let as_ = a * a; + let mut vi = rrtpi * a * (-0.5 * ah * ah).exp(); + let mut zi = znorm1(ah) / h; + let y = 1.0 / hs; + loop { + tf += zi * c2[i - 1]; + if i > m { + tf *= rrtpi * (-0.5 * hs).exp(); + return tf; + } + zi = y * ((ii as f64) * zi - vi); + vi *= as_; + i += 1; + ii += 2; } - zi = y * ((ii as f64) * zi - vi); - vi *= as_; - i += 1; - ii += 2; } - } - 4 => { - let maxii = m + m + 1; - let mut ii = 1; - let mut tf = 0.0; - let hs = h * h; - let as_ = -a * a; - let mut ai = rtwopi * a * (-0.5 * hs * (1.0 - as_)).exp(); - let mut yi = 1.0; - loop { - tf += ai * yi; - if ii >= maxii { - return tf; + 4 => { + let maxii = m + m + 1; + let mut ii = 1; + let mut tf = 0.0; + let hs = h * h; + let as_ = -a * a; + let mut ai = rtwopi * a * (-0.5 * hs * (1.0 - as_)).exp(); + let mut yi = 1.0; + loop { + tf += ai * yi; + if ii >= maxii { + return tf; + } + ii += 2; + yi = (1.0 - hs * yi) / (ii as f64); + ai *= as_; } - ii += 2; - yi = (1.0 - hs * yi) / (ii as f64); - ai *= as_; } - } - 5 => { - let mut tf = 0.0; - let as_ = a * a; - let hs = -0.5 * h * h; - for i in 0..m { - let r = 1.0 + as_ * pts[i]; - tf += wts[i] * (hs * r).exp() / r; + 5 => { + let mut tf = 0.0; + let as_ = a * a; + let hs = -0.5 * h * h; + for i in 0..m { + let r = 1.0 + as_ * pts[i]; + tf += wts[i] * (hs * r).exp() / r; + } + tf *= a; + tf } - tf *= a; - tf - } - 6 => { - let normh = znorm2(h); - let mut tf = 0.5 * normh * (1.0 - normh); - let y = 1.0 - a; - let r = (y / (1.0 + a)).atan(); - if r != 0.0 { - tf -= rtwopi * r * (-0.5 * y * h * h / r).exp(); + 6 => { + let normh = znorm2(h); + let mut tf = 0.5 * normh * (1.0 - normh); + let y = 1.0 - a; + let r = (y / (1.0 + a)).atan(); + if r != 0.0 { + tf -= rtwopi * r * (-0.5 * y * h * h / r).exp(); + } + tf } - tf + _ => 0.0, } - _ => 0.0, } -} -fn phi(x: f64) -> f64 { - normal_cdf(x, 0.0, 1.0) -} + // P(0 ≤ Z ≤ x) + fn znorm1(x: f64) -> f64 { + phi(x) - 0.5 + } + + // P(x ≤ Z < ∞) + fn znorm2(x: f64) -> f64 { + 1.0 - phi(x) + } -// P(0 ≤ Z ≤ x) -fn znorm1(x: f64) -> f64 { - phi(x) - 0.5 + t } -// P(x ≤ Z < ∞) -fn znorm2(x: f64) -> f64 { - 1.0 - phi(x) +/// standard normal cdf +fn phi(x: f64) -> f64 { + normal_cdf(x, 0.0, 1.0) } From 7ca4e29de4c520a1287316b39f10b0a0b7f1868e Mon Sep 17 00:00:00 2001 From: JamboChen Date: Fri, 18 Oct 2024 00:54:03 +0800 Subject: [PATCH 23/24] Collect MIT License functions --- rand_distr/tests/cdf.rs | 115 +---------------------- rand_distr/tests/common/mod.rs | 1 + rand_distr/tests/common/spec_func.rs | 134 +++++++++++++++++++++++++++ 3 files changed, 139 insertions(+), 111 deletions(-) create mode 100644 rand_distr/tests/common/mod.rs create mode 100644 rand_distr/tests/common/spec_func.rs diff --git a/rand_distr/tests/cdf.rs b/rand_distr/tests/cdf.rs index fde345a124..8eb22740e2 100644 --- a/rand_distr/tests/cdf.rs +++ b/rand_distr/tests/cdf.rs @@ -6,6 +6,7 @@ // option. This file may not be copied, modified, or distributed // except according to those terms. +mod common; use core::f64; use num_traits::AsPrimitive; @@ -353,6 +354,7 @@ fn frechet() { #[test] fn zeta() { fn cdf(k: i64, s: f64) -> f64 { + use common::spec_func::zeta_func; if k < 1 { return 0.0; } @@ -609,6 +611,8 @@ fn hypergeometric() { fn poisson() { use rand_distr::Poisson; fn cdf(k: i64, lambda: f64) -> f64 { + use common::spec_func::gamma_lr; + if k < 0 || lambda <= 0.0 { return 0.0; } @@ -635,79 +639,6 @@ fn ln_binomial(n: u64, k: u64) -> f64 { ln_factorial(n) - ln_factorial(k) - ln_factorial(n - k) } -/// https://docs.rs/statrs/latest/src/statrs/function/gamma.rs.html#260-347 -fn gamma_lr(a: f64, x: f64) -> f64 { - let eps = 0.000000000000001; - let big = 4503599627370496.0; - let big_inv = 2.220_446_049_250_313e-16; - - let ax = a * x.ln() - x - a.lgamma().0; - if ax < -709.782_712_893_384 { - if a < x { - return 1.0; - } - return 0.0; - } - if x <= 1.0 || x <= a { - let mut r2 = a; - let mut c2 = 1.0; - let mut ans2 = 1.0; - loop { - r2 += 1.0; - c2 *= x / r2; - ans2 += c2; - - if c2 / ans2 <= eps { - break; - } - } - return ax.exp() * ans2 / a; - } - - let mut y = 1.0 - a; - let mut z = x + y + 1.0; - let mut c = 0; - - let mut p3 = 1.0; - let mut q3 = x; - let mut p2 = x + 1.0; - let mut q2 = z * x; - let mut ans = p2 / q2; - - loop { - y += 1.0; - z += 2.0; - c += 1; - let yc = y * f64::from(c); - - let p = p2 * z - p3 * yc; - let q = q2 * z - q3 * yc; - - p3 = p2; - p2 = p; - q3 = q2; - q2 = q; - - if p.abs() > big { - p3 *= big_inv; - p2 *= big_inv; - q3 *= big_inv; - q2 *= big_inv; - } - - if q != 0.0 { - let nextans = p / q; - let error = ((ans - nextans) / nextans).abs(); - ans = nextans; - - if error <= eps { - break; - } - } - } - 1.0 - ax.exp() * ans -} - fn gen_harmonic(n: u64, m: f64) -> f64 { match n { 0 => 1.0, @@ -715,44 +646,6 @@ fn gen_harmonic(n: u64, m: f64) -> f64 { } } -/// https://docs.rs/spfunc/latest/src/spfunc/zeta.rs.html#20-43 -fn zeta_func(s: f64) -> f64 { - if s == 1.0 { - return f64::INFINITY; - } - - let mut akn = vec![1.0]; - let mut two_pow = 0.5; - let head = 1.0 / (1.0 - 2.0.powf(1.0 - s)); - let mut tail_prev = 0.0; - let mut tail = two_pow * akn[0]; - - while (tail - tail_prev).abs() >= f64::EPSILON { - update_akn(&mut akn, s); - two_pow /= 2.0; - tail_prev = tail; - tail += two_pow * akn.iter().sum::(); - } - - head * tail -} - -fn update_akn(akn: &mut Vec, s: f64) { - let n = akn.len() - 1; - - let n1 = n as f64 + 1.0; - - akn.iter_mut().enumerate().for_each(|(k, a)| { - let num = n1; - let den = n + 1 - k; - *a *= num / den as f64; - }); - let p1 = -1.0 / n1; - let p2 = (n1 / (n1 + 1.0)).powf(s); - - akn.push(p1 * p2 * akn[n]); -} - /// [1] Patefield, M. (2000). Fast and Accurate Calculation of Owen’s T Function. /// Journal of Statistical Software, 5(5), 1–25. /// https://doi.org/10.18637/jss.v005.i05 diff --git a/rand_distr/tests/common/mod.rs b/rand_distr/tests/common/mod.rs new file mode 100644 index 0000000000..4fa8db7de6 --- /dev/null +++ b/rand_distr/tests/common/mod.rs @@ -0,0 +1 @@ +pub mod spec_func; diff --git a/rand_distr/tests/common/spec_func.rs b/rand_distr/tests/common/spec_func.rs new file mode 100644 index 0000000000..6515761082 --- /dev/null +++ b/rand_distr/tests/common/spec_func.rs @@ -0,0 +1,134 @@ +// MIT License + +// Copyright (c) 2016 Michael Ma + +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: + +// The above copyright notice and this permission notice shall be included in all +// copies or substantial portions of the Software. + +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. + +/// https://docs.rs/statrs/latest/src/statrs/function/gamma.rs.html#260-347 +use special::Primitive; + +pub fn gamma_lr(a: f64, x: f64) -> f64 { + let eps = 0.000000000000001; + let big = 4503599627370496.0; + let big_inv = 2.220_446_049_250_313e-16; + + let ax = a * x.ln() - x - a.lgamma().0; + if ax < -709.782_712_893_384 { + if a < x { + return 1.0; + } + return 0.0; + } + if x <= 1.0 || x <= a { + let mut r2 = a; + let mut c2 = 1.0; + let mut ans2 = 1.0; + loop { + r2 += 1.0; + c2 *= x / r2; + ans2 += c2; + + if c2 / ans2 <= eps { + break; + } + } + return ax.exp() * ans2 / a; + } + + let mut y = 1.0 - a; + let mut z = x + y + 1.0; + let mut c = 0; + + let mut p3 = 1.0; + let mut q3 = x; + let mut p2 = x + 1.0; + let mut q2 = z * x; + let mut ans = p2 / q2; + + loop { + y += 1.0; + z += 2.0; + c += 1; + let yc = y * f64::from(c); + + let p = p2 * z - p3 * yc; + let q = q2 * z - q3 * yc; + + p3 = p2; + p2 = p; + q3 = q2; + q2 = q; + + if p.abs() > big { + p3 *= big_inv; + p2 *= big_inv; + q3 *= big_inv; + q2 *= big_inv; + } + + if q != 0.0 { + let nextans = p / q; + let error = ((ans - nextans) / nextans).abs(); + ans = nextans; + + if error <= eps { + break; + } + } + } + 1.0 - ax.exp() * ans +} + +/// https://docs.rs/spfunc/latest/src/spfunc/zeta.rs.html#20-43 +pub fn zeta_func(s: f64) -> f64 { + fn update_akn(akn: &mut Vec, s: f64) { + let n = akn.len() - 1; + + let n1 = n as f64 + 1.0; + + akn.iter_mut().enumerate().for_each(|(k, a)| { + let num = n1; + let den = n + 1 - k; + *a *= num / den as f64; + }); + let p1 = -1.0 / n1; + let p2 = (n1 / (n1 + 1.0)).powf(s); + + akn.push(p1 * p2 * akn[n]); + } + + if s == 1.0 { + return f64::INFINITY; + } + + let mut akn = vec![1.0]; + let mut two_pow = 0.5; + let head = 1.0 / (1.0 - 2.0.powf(1.0 - s)); + let mut tail_prev = 0.0; + let mut tail = two_pow * akn[0]; + + while (tail - tail_prev).abs() >= f64::EPSILON { + update_akn(&mut akn, s); + two_pow /= 2.0; + tail_prev = tail; + tail += two_pow * akn.iter().sum::(); + } + + head * tail +} From b359a588a34a0d7a7aa17e9049313b41feedd267 Mon Sep 17 00:00:00 2001 From: JamboChen Date: Fri, 18 Oct 2024 01:03:58 +0800 Subject: [PATCH 24/24] Add `assert` for `gamma_lr` --- rand_distr/tests/common/spec_func.rs | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/rand_distr/tests/common/spec_func.rs b/rand_distr/tests/common/spec_func.rs index 6515761082..bd39a79fb4 100644 --- a/rand_distr/tests/common/spec_func.rs +++ b/rand_distr/tests/common/spec_func.rs @@ -24,6 +24,27 @@ use special::Primitive; pub fn gamma_lr(a: f64, x: f64) -> f64 { + if a.is_nan() || x.is_nan() { + return f64::NAN; + } + + if a <= 0.0 || a == f64::INFINITY { + panic!("a must be positive and finite"); + } + if x <= 0.0 || x == f64::INFINITY { + panic!("x must be positive and finite"); + } + + const ACC: f64 = 0.0000000000000011102230246251565; + + if a.abs() < ACC { + return 1.0; + } + + if x.abs() < ACC { + return 0.0; + } + let eps = 0.000000000000001; let big = 4503599627370496.0; let big_inv = 2.220_446_049_250_313e-16;