Fix catastrophic cancellation error in poisson logpdf when k and mu are large

[rahulb0802]

This replaces the current direct-formula method in logpdf with Loader's (2000) saddle-point approx, ported from a change I've been working on in SciPy's scipy.stats.poisson.logpmf. The current formula loses precision due to cancellation when k is close to mu and both are large. The new form implemented avoids this.

Two new helpers are added to poisson_detail:

• stirlerr(n): Stirling's series remainder. Direct path for small n, asymptotic series for large n
• bd0(mean, k): Deviance term. Direct path and series expansion for |v| < 0.1 to avoid cancellation when k and mu are close

All existing poisson tests pass (across all real types). I also added accuracy checks for other values, validated against mpmath (1000-digit) reference values, including large-magnitude cases. No regressions as per CI checks.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 95.40%. Comparing base (64a8d75) to head (487009b).

Additional details and impacted files

@@           Coverage Diff            @@
##           develop    #1438   +/-   ##
========================================
  Coverage    95.39%   95.40%           
========================================
  Files          829      829           
  Lines        69186    69219   +33     
========================================
+ Hits         66003    66039   +36     
+ Misses        3183     3180    -3     

Files with missing lines	Coverage Δ
include/boost/math/distributions/poisson.hpp	86.03% <100.00%> (+2.16%)	⬆️
test/test_poisson.cpp	99.19% <100.00%> (+0.06%)	⬆️

... and 1 file with indirect coverage changes

---

Continue to review full report in Codecov by Harness.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 64a8d75...487009b. Read the comment docs.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[jzmaddock]

Thanks for this @rahulb0802 , this is a good idea, but not quite right just yet. The problem is that the new routines evaluate a fixed number of terms of the Stirling remainder, and therefore are definitely not arbitrary precision. BTW your new test values are limited to double precision too, so the tests should NOT be passing for any long double larger than 64 bits, certainly not for the 128-bit ones, I need to investigate why that is.

But for now taking the log of the PDF as the true result, and with old_logpdf a copy of the previous logpdf code:

    using boost::math::poisson_distribution;

    using mp = boost::multiprecision::number<boost::multiprecision::debug_adaptor<boost::multiprecision::backends::cpp_bin_float<200>>>;

    std::int64_t mean = 20;
    std::int64_t k    =  18;

    std::cout << std::setprecision(18);

    mp val = log(pdf(poisson_distribution<mp>(mean), mp(k)));
    double dval = (double)val;

    std::cout << "mean =             " << mean << " k = " << k << std::endl;
    std::cout << "MP logpdf:         " << boost::math::relative_difference(logpdf(poisson_distribution<mp>(mean), mp(k)), val) << std::endl;
    std::cout << "MP old logpdf:     " << boost::math::relative_difference(old_logpdf(poisson_distribution<mp>(mean), mp(k)), val) << std::endl;
    std::cout << "double logpdf:     " << boost::math::relative_difference(logpdf(poisson_distribution<double>(mean), double(k)), dval) << std::endl;
    std::cout << "double old logpdf: " << boost::math::relative_difference(old_logpdf(poisson_distribution<double>(mean), double(k)), dval) << std::endl;

Gives:

mean =             20 k = 18
MP logpdf:         1.19456007034065095e-17
MP old logpdf:     1.03043973124271719e-198
ddouble logpdf:     1.79628534341220374e-16
double old logpdf: 8.98142671706101575e-16

So when evaluated at 200-decimal precision the new routine gives just 17 decimal places compared to nearly full precision before. Evaluating at larger mean looks a little better:

mean =             1000000 k = 1300000
MP logpdf:         2.60446354798830633e-75
MP old logpdf:     1.30516145170170221e-197
double logpdf:     2.12532375379900552e-15
double old logpdf: 5.75608516653897354e-14

But we're further out into the asymptotic region for Sterlings formula so fewer terms generate more digits of precision.

So... what to do? I can't see Sterling converging to full arbitrary precision for small means (from memory the terms become alternating and divergent eventually?) Although we might get there if the is_small check is appropriately tweaked inside stirlerr? We do have an arbitrary precision version of Sterling here:

math/include/boost/math/special_functions/gamma.hpp

Line 490 in 9b192b0

T scaled_tgamma_no_lanczos(const T& z, const Policy& pol, bool islog = false)

which calculates tgamma(z) / (z/e)^z, could that be pressed into action here?

[rahulb0802]

Yep, that is the fix, thanks for pointing this out. The sum variable being accumulated in this function is literally the exact same thing I was trying to compute in stirlerr, but not capped at a fixed number of terms (five in my original implementation). So I moved that loop part to its own function and now the scaled_tgamma_no_lanczos function just calls that before performing final computations.

bd0 had a similar problem but the fix was much simpler there. Instead of using 10 terms, just put in a convergence check. But doing it this way introduces an edge case for nonconvergence when k = mean. However, since we know the series is trivially zero at that point, we can safely skip the loop.

Running your original case:

using mp = boost::multiprecision::number<boost::multiprecision::backends::cpp_bin_float<200>>;

std::cout << std::setprecision(18);

std::int64_t mean = 20;
std::int64_t k    = 18;

mp val = log(boost::math::pdf(poisson_distribution<mp>(mean), mp(k)));
double dval = (double)val;

std::cout << "mean =             " << mean << " k = " << k << std::endl;
std::cout << "MP logpdf:         " << boost::math::relative_difference(boost::math::logpdf(poisson_distribution<mp>(mean), mp(k)), val) << std::endl;
std::cout << "double logpdf:     " << boost::math::relative_difference(boost::math::logpdf(poisson_distribution<double>(mean), double(k)), dval) << std::endl;

mean =             20 k = 18
MP logpdf:         6.86959820828478126e-200
double logpdf:     1.79628534341220374e-16

mean =             1000000 k = 1300000
MP logpdf:         8.80528085179990715e-199
double logpdf:     2.12532375379900552e-15

Double results remain unchanged but the mp result is substantially better in terms of precision.