From 9cdc9847a64a5424f695c241cf38c9419852b775 Mon Sep 17 00:00:00 2001 From: rahulb0802 Date: Tue, 25 Aug 2026 21:49:29 -0500 Subject: [PATCH 1/3] Ports Loader (2000) SP approx to prevent cancellation in extreme parameter regimes --- include/boost/math/distributions/poisson.hpp | 48 ++++++++++++++- test/test_poisson.cpp | 61 +++++++++++++++++++- 2 files changed, 107 insertions(+), 2 deletions(-) diff --git a/include/boost/math/distributions/poisson.hpp b/include/boost/math/distributions/poisson.hpp index c5e4404335..e8a31c6b62 100644 --- a/include/boost/math/distributions/poisson.hpp +++ b/include/boost/math/distributions/poisson.hpp @@ -50,6 +50,7 @@ #include // factorials. #include // for root finding. #include +#include namespace boost { @@ -141,6 +142,51 @@ namespace boost return true; } // bool check_dist_and_prob + template + BOOST_MATH_GPU_ENABLED inline RealType stirlerr(const RealType& n) { + BOOST_MATH_STD_USING // for ADL of std functions. + using boost::math::lgamma; + + const RealType S0 = RealType(1)/12; + const RealType S1 = RealType(1)/360; + const RealType S2 = RealType(1)/1260; + const RealType S3 = RealType(1)/1680; + const RealType S4 = RealType(1)/1188; + + bool is_small = n < 15; + if (is_small) { + return lgamma(n + 1) - (n * log(n) - n + 0.5 * log(2 * boost::math::constants::pi() * n)); + } else { + RealType n2 = n * n; + return (S0 - (S1 - (S2 - (S3 - S4/n2)/n2)/n2)/n2)/n; + } + + } + + template + BOOST_MATH_GPU_ENABLED inline RealType bd0(const RealType& mean, const RealType& k) { + BOOST_MATH_STD_USING // for ADL of std functions. + + bool is_close = abs(k - mean) < RealType(0.1) * (k + mean); + + if (is_close) { + RealType v = (k - mean) / (k + mean); + RealType v2 = v * v; + RealType series_term = ((k - mean) * (k - mean)) / (k + mean); + + RealType term = 2 * k * v; + for (int i = 1; i < 11; ++i) { + term *= v2; + series_term += term / (2 * i + 1); + } + return series_term; + } else { + RealType direct = (k == 0) ? RealType(0) : k * log(k / mean) + mean - k; + return direct; + } + + } + } // namespace poisson_detail BOOST_MATH_EXPORT template > @@ -304,7 +350,7 @@ namespace boost // Special case where k and lambda are both positive if(k > 0 && mean > 0) { - return -lgamma(k+1) + k*log(mean) - mean; + return -poisson_detail::stirlerr(k) - poisson_detail::bd0(mean, k) - RealType(0.5) * log(2 * boost::math::constants::pi() * k); } result = log(pdf(dist, k)); diff --git a/test/test_poisson.cpp b/test/test_poisson.cpp index 96e5f12d73..8a55028f21 100644 --- a/test/test_poisson.cpp +++ b/test/test_poisson.cpp @@ -244,7 +244,66 @@ void test_spots(RealType) static_cast(20)), // K>> mean log(static_cast(8.277463646553730E-009)), // probability. tolerance); - + + BOOST_CHECK_CLOSE( + logpdf(poisson_distribution(static_cast(14)), // mean 14. + static_cast(14)), + static_cast(-2.244418568125061), // probability. + tolerance); + + BOOST_CHECK_CLOSE( + logpdf(poisson_distribution(static_cast(20)), // mean 20. + static_cast(18)), + static_cast(-2.472264284061216), // probability. + tolerance); + + // Cases below require around 15+ significant decimal digits to represent + // k / mean meaningfully, so skip for float. + if (std::numeric_limits::digits10 > 15) + { + BOOST_CHECK_CLOSE( + logpdf(poisson_distribution(static_cast(1000000)), + static_cast(1300000)), + static_cast(-41081.501683746894), + tolerance); + + BOOST_CHECK_CLOSE( + logpdf(poisson_distribution(static_cast(1e8)), + static_cast(8e7)), + static_cast(-2148525.91257035), + tolerance); + + BOOST_CHECK_CLOSE( + logpdf(poisson_distribution(static_cast(1e10)), + static_cast(105e9)), + static_cast(-151894402015.7727), + tolerance); + + BOOST_CHECK_CLOSE( + logpdf(poisson_distribution(static_cast(1e15)), + static_cast(8e14)), // |v| > 0.1 boundary + static_cast(-21485158948650.273), + tolerance); + + BOOST_CHECK_CLOSE( + logpdf(poisson_distribution(static_cast(1e15)), + static_cast(12e14)), // |v| < 0.1 boundary + static_cast(-18785868152763.832), + tolerance); + + BOOST_CHECK_CLOSE( + logpdf(poisson_distribution(static_cast(1e16)), + static_cast(1e16)), // old formula returns 0.0 here + static_cast(-19.339619277157038), + tolerance); + + BOOST_CHECK_CLOSE( + logpdf(poisson_distribution(static_cast(5e15)), + static_cast(5e15)), + static_cast(-18.993045686877064), + tolerance); + } + // CDF BOOST_CHECK_CLOSE( cdf(poisson_distribution(static_cast(1)), // mean unity. From 82a1da8785b4986be36dba00f0d07bd07d01d7c5 Mon Sep 17 00:00:00 2001 From: rahulb0802 Date: Tue, 25 Aug 2026 21:55:36 -0500 Subject: [PATCH 2/3] Add comments to approximation code and testing --- include/boost/math/distributions/poisson.hpp | 7 ++++++- test/test_poisson.cpp | 18 +++++++++--------- 2 files changed, 15 insertions(+), 10 deletions(-) diff --git a/include/boost/math/distributions/poisson.hpp b/include/boost/math/distributions/poisson.hpp index e8a31c6b62..68c0463eb6 100644 --- a/include/boost/math/distributions/poisson.hpp +++ b/include/boost/math/distributions/poisson.hpp @@ -147,12 +147,14 @@ namespace boost BOOST_MATH_STD_USING // for ADL of std functions. using boost::math::lgamma; + // Stirling's series coefficients const RealType S0 = RealType(1)/12; const RealType S1 = RealType(1)/360; const RealType S2 = RealType(1)/1260; const RealType S3 = RealType(1)/1680; const RealType S4 = RealType(1)/1188; + // Use Stirling's series if n is small; use the direct formula otherwise bool is_small = n < 15; if (is_small) { return lgamma(n + 1) - (n * log(n) - n + 0.5 * log(2 * boost::math::constants::pi() * n)); @@ -167,9 +169,10 @@ namespace boost BOOST_MATH_GPU_ENABLED inline RealType bd0(const RealType& mean, const RealType& k) { BOOST_MATH_STD_USING // for ADL of std functions. + // Calculate v = (k - mean) / (k + mean) from Loader (2000) approximation bool is_close = abs(k - mean) < RealType(0.1) * (k + mean); - if (is_close) { + if (is_close) { // Use the series approximation if |v| < 0.1 RealType v = (k - mean) / (k + mean); RealType v2 = v * v; RealType series_term = ((k - mean) * (k - mean)) / (k + mean); @@ -181,6 +184,7 @@ namespace boost } return series_term; } else { + // Use the direct formula if |v| >= 0.1 RealType direct = (k == 0) ? RealType(0) : k * log(k / mean) + mean - k; return direct; } @@ -350,6 +354,7 @@ namespace boost // Special case where k and lambda are both positive if(k > 0 && mean > 0) { + // Use the Loader (2000) saddle-point approximation for logpdf calculation return -poisson_detail::stirlerr(k) - poisson_detail::bd0(mean, k) - RealType(0.5) * log(2 * boost::math::constants::pi() * k); } diff --git a/test/test_poisson.cpp b/test/test_poisson.cpp index 8a55028f21..9147e6787f 100644 --- a/test/test_poisson.cpp +++ b/test/test_poisson.cpp @@ -248,13 +248,13 @@ void test_spots(RealType) BOOST_CHECK_CLOSE( logpdf(poisson_distribution(static_cast(14)), // mean 14. static_cast(14)), - static_cast(-2.244418568125061), // probability. + static_cast(-2.244418568125061), // probability (already in log space). tolerance); BOOST_CHECK_CLOSE( logpdf(poisson_distribution(static_cast(20)), // mean 20. static_cast(18)), - static_cast(-2.472264284061216), // probability. + static_cast(-2.472264284061216), // probability (already in log space). tolerance); // Cases below require around 15+ significant decimal digits to represent @@ -262,43 +262,43 @@ void test_spots(RealType) if (std::numeric_limits::digits10 > 15) { BOOST_CHECK_CLOSE( - logpdf(poisson_distribution(static_cast(1000000)), + logpdf(poisson_distribution(static_cast(1000000)), // mean 1000000. static_cast(1300000)), static_cast(-41081.501683746894), tolerance); BOOST_CHECK_CLOSE( - logpdf(poisson_distribution(static_cast(1e8)), + logpdf(poisson_distribution(static_cast(1e8)), // mean 1e8. static_cast(8e7)), static_cast(-2148525.91257035), tolerance); BOOST_CHECK_CLOSE( - logpdf(poisson_distribution(static_cast(1e10)), + logpdf(poisson_distribution(static_cast(1e10)), // mean 1e10. static_cast(105e9)), static_cast(-151894402015.7727), tolerance); BOOST_CHECK_CLOSE( - logpdf(poisson_distribution(static_cast(1e15)), + logpdf(poisson_distribution(static_cast(1e15)), // mean 1e15. static_cast(8e14)), // |v| > 0.1 boundary static_cast(-21485158948650.273), tolerance); BOOST_CHECK_CLOSE( - logpdf(poisson_distribution(static_cast(1e15)), + logpdf(poisson_distribution(static_cast(1e15)), // mean 1e15. static_cast(12e14)), // |v| < 0.1 boundary static_cast(-18785868152763.832), tolerance); BOOST_CHECK_CLOSE( - logpdf(poisson_distribution(static_cast(1e16)), + logpdf(poisson_distribution(static_cast(1e16)), // mean 1e16. static_cast(1e16)), // old formula returns 0.0 here static_cast(-19.339619277157038), tolerance); BOOST_CHECK_CLOSE( - logpdf(poisson_distribution(static_cast(5e15)), + logpdf(poisson_distribution(static_cast(5e15)), // mean 5e15. static_cast(5e15)), static_cast(-18.993045686877064), tolerance); From 487009bc32d72892c68fe1331905f2351973cde7 Mon Sep 17 00:00:00 2001 From: rahulb0802 Date: Tue, 25 Aug 2026 22:06:51 -0500 Subject: [PATCH 3/3] mpmath reference in testing comment --- test/test_poisson.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/test_poisson.cpp b/test/test_poisson.cpp index 9147e6787f..22daed882f 100644 --- a/test/test_poisson.cpp +++ b/test/test_poisson.cpp @@ -245,6 +245,8 @@ void test_spots(RealType) log(static_cast(8.277463646553730E-009)), // probability. tolerance); + // New test cases for Loader (2000) saddle-point approximation. Probs already + // in log space. Values calculated using mpmath (1000-digit precision). BOOST_CHECK_CLOSE( logpdf(poisson_distribution(static_cast(14)), // mean 14. static_cast(14)),