From 24786e6468dd3bc340926781f9366b3d8d710956 Mon Sep 17 00:00:00 2001 From: Benjamin Lieser Date: Sat, 12 Oct 2024 15:34:33 +0200 Subject: [PATCH 1/9] Fix bug in hypergeometric rejection sampling --- rand_distr/src/hypergeometric.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rand_distr/src/hypergeometric.rs b/rand_distr/src/hypergeometric.rs index c1f1d4ef234..4f5ac3429b6 100644 --- a/rand_distr/src/hypergeometric.rs +++ b/rand_distr/src/hypergeometric.rs @@ -359,7 +359,7 @@ impl Distribution for Hypergeometric { } else { for i in (y as u64 + 1)..=(m as u64) { f *= i as f64 * (n2 - k + i) as f64; - f /= (n1 - i) as f64 * (k - i) as f64; + f /= (n1 - i + 1) as f64 * (k - i + 1) as f64; } } From 3ee2a3c2b1fe1784e601081bdeabf82be4ad47b5 Mon Sep 17 00:00:00 2001 From: Benjamin Lieser Date: Sat, 12 Oct 2024 15:39:27 +0200 Subject: [PATCH 2/9] better stirling --- rand_distr/src/hypergeometric.rs | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/rand_distr/src/hypergeometric.rs b/rand_distr/src/hypergeometric.rs index 4f5ac3429b6..fad878275ae 100644 --- a/rand_distr/src/hypergeometric.rs +++ b/rand_distr/src/hypergeometric.rs @@ -131,10 +131,12 @@ fn fraction_of_products_of_factorials(numerator: (u64, u64), denominator: (u64, result } +const LOGSQRT2PI: f64 = 0.91893853320467274178; // log(sqrt(2)*pi) + fn ln_of_factorial(v: f64) -> f64 { // the paper calls for ln(v!), but also wants to pass in fractions, // so we need to use Stirling's approximation to fill in the gaps: - v * v.ln() - v + (v + 0.5) * v.ln() - v + LOGSQRT2PI + 1.0 / (12.0 * v) } impl Hypergeometric { From 76900b1f77e3835ed7450c8d49d4317370c91b67 Mon Sep 17 00:00:00 2001 From: Benjamin Lieser Date: Sat, 12 Oct 2024 15:58:53 +0200 Subject: [PATCH 3/9] more precise stirling with test --- rand_distr/src/hypergeometric.rs | 20 +++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/rand_distr/src/hypergeometric.rs b/rand_distr/src/hypergeometric.rs index fad878275ae..d18202c96f1 100644 --- a/rand_distr/src/hypergeometric.rs +++ b/rand_distr/src/hypergeometric.rs @@ -136,7 +136,12 @@ const LOGSQRT2PI: f64 = 0.91893853320467274178; // log(sqrt(2)*pi) fn ln_of_factorial(v: f64) -> f64 { // the paper calls for ln(v!), but also wants to pass in fractions, // so we need to use Stirling's approximation to fill in the gaps: - (v + 0.5) * v.ln() - v + LOGSQRT2PI + 1.0 / (12.0 * v) + + // shift v by 3, because Stirling is bad for small values + let v_3 = v + 3.0; + let ln_fac = (v_3 + 0.5) * v_3.ln() - v_3 + LOGSQRT2PI + 1.0 / (12.0 * v_3); + // make the correction for the shift + ln_fac - ((v + 3.0) * (v + 2.0) * (v + 1.0)).ln() } impl Hypergeometric { @@ -443,6 +448,8 @@ impl Distribution for Hypergeometric { #[cfg(test)] mod test { + use std::{dbg, vec}; + use super::*; #[test] @@ -496,4 +503,15 @@ mod test { fn hypergeometric_distributions_can_be_compared() { assert_eq!(Hypergeometric::new(1, 2, 3), Hypergeometric::new(1, 2, 3)); } + + #[test] + fn stirling() { + use special::Gamma; + let test = vec![0.5, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]; + for &v in test.iter() { + let ln_fac = ln_of_factorial(v); + dbg!(ln_fac); + assert!((ln_fac - (v + 1.0).ln_gamma().0).abs() < 1e-4); + } + } } From 23d26ddd3b9e10756744f6ee7608b632740855a4 Mon Sep 17 00:00:00 2001 From: Benjamin Lieser Date: Sat, 12 Oct 2024 16:06:35 +0200 Subject: [PATCH 4/9] clippy --- rand_distr/src/hypergeometric.rs | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/rand_distr/src/hypergeometric.rs b/rand_distr/src/hypergeometric.rs index d18202c96f1..6f7bc4d0af3 100644 --- a/rand_distr/src/hypergeometric.rs +++ b/rand_distr/src/hypergeometric.rs @@ -506,12 +506,11 @@ mod test { #[test] fn stirling() { - use special::Gamma; - let test = vec![0.5, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]; + let test = [0.5, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]; for &v in test.iter() { let ln_fac = ln_of_factorial(v); dbg!(ln_fac); - assert!((ln_fac - (v + 1.0).ln_gamma().0).abs() < 1e-4); + assert!((special::Gamma::ln_gamma(v + 1.0).0 - ln_fac).abs() < 1e-4); } } } From 7644a1c6b71236e4bd99ebbd7e78a3ca2597cb7b Mon Sep 17 00:00:00 2001 From: Benjamin Lieser Date: Sat, 12 Oct 2024 16:09:31 +0200 Subject: [PATCH 5/9] remove dbg --- rand_distr/src/hypergeometric.rs | 2 -- 1 file changed, 2 deletions(-) diff --git a/rand_distr/src/hypergeometric.rs b/rand_distr/src/hypergeometric.rs index 6f7bc4d0af3..71a5855a98a 100644 --- a/rand_distr/src/hypergeometric.rs +++ b/rand_distr/src/hypergeometric.rs @@ -448,7 +448,6 @@ impl Distribution for Hypergeometric { #[cfg(test)] mod test { - use std::{dbg, vec}; use super::*; @@ -509,7 +508,6 @@ mod test { let test = [0.5, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]; for &v in test.iter() { let ln_fac = ln_of_factorial(v); - dbg!(ln_fac); assert!((special::Gamma::ln_gamma(v + 1.0).0 - ln_fac).abs() < 1e-4); } } From 35b72e6535cd26bcf4decbda011e4d6f29681c4c Mon Sep 17 00:00:00 2001 From: Benjamin Lieser Date: Sat, 12 Oct 2024 16:13:07 +0200 Subject: [PATCH 6/9] value stability --- rand_distr/tests/value_stability.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rand_distr/tests/value_stability.rs b/rand_distr/tests/value_stability.rs index b142741e77c..330119b68f6 100644 --- a/rand_distr/tests/value_stability.rs +++ b/rand_distr/tests/value_stability.rs @@ -105,7 +105,7 @@ fn hypergeometric_stability() { test_samples( 7221, Hypergeometric::new(100, 50, 50).unwrap(), - &[23, 27, 26, 27, 22, 24, 31, 22], + &[23, 27, 26, 27, 22, 25, 31, 25], ); // Algorithm H2PE } From 7ef4df30eafcc226c3b8ff91095511258297d7df Mon Sep 17 00:00:00 2001 From: Benjamin Lieser Date: Thu, 24 Oct 2024 13:48:09 +0200 Subject: [PATCH 7/9] changelog --- rand_distr/CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/rand_distr/CHANGELOG.md b/rand_distr/CHANGELOG.md index c9ab729be8c..87e2f68f1f0 100644 --- a/rand_distr/CHANGELOG.md +++ b/rand_distr/CHANGELOG.md @@ -10,6 +10,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - 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) +- Fix bug in `Hypergeometric`, this is a Value-breaking change (#1510) ### Added - Add plots for `rand_distr` distributions to documentation (#1434) From a9e3945fef2755b707da7b754684908ca1f26a73 Mon Sep 17 00:00:00 2001 From: Benjamin Lieser Date: Thu, 24 Oct 2024 13:55:36 +0200 Subject: [PATCH 8/9] add prev failing cdf test --- rand_distr/tests/cdf.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rand_distr/tests/cdf.rs b/rand_distr/tests/cdf.rs index 8eb22740e2b..10d04e21563 100644 --- a/rand_distr/tests/cdf.rs +++ b/rand_distr/tests/cdf.rs @@ -598,7 +598,7 @@ fn hypergeometric() { (60, 10, 7), (70, 20, 50), (100, 50, 10), - // (100, 50, 49), // Fail case + (100, 50, 49), ]; for (seed, (n, k, n_)) in parameters.into_iter().enumerate() { From d96e2195274eb9451702d1682c016605842e2f8f Mon Sep 17 00:00:00 2001 From: Benjamin Lieser Date: Wed, 13 Nov 2024 11:53:57 +0100 Subject: [PATCH 9/9] fix comment --- rand_distr/src/hypergeometric.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rand_distr/src/hypergeometric.rs b/rand_distr/src/hypergeometric.rs index 71a5855a98a..f446357530b 100644 --- a/rand_distr/src/hypergeometric.rs +++ b/rand_distr/src/hypergeometric.rs @@ -131,7 +131,7 @@ fn fraction_of_products_of_factorials(numerator: (u64, u64), denominator: (u64, result } -const LOGSQRT2PI: f64 = 0.91893853320467274178; // log(sqrt(2)*pi) +const LOGSQRT2PI: f64 = 0.91893853320467274178; // log(sqrt(2*pi)) fn ln_of_factorial(v: f64) -> f64 { // the paper calls for ln(v!), but also wants to pass in fractions,