Repository navigation
BUG: Variance inflation in random_poisson (and dependent samplers) for large lam (≳1e15) due to float64 rounding #31986
Description
Describe the issue:
np.random.default_rng(0).poisson(lam) (as well as the legacy API
np.random.poisson) draws samples whose empirical variance is substantially
larger than the theoretical variance (which should be lam), once lam is large.
The sample mean stays correct, but var(x) / lam grows well above 1.0.
Measured var/lam on current main / 2.4.4 (size = 2_000_000):
| lam | var / lam |
|---|---|
| 1e10 | ~1.000 |
| 1e14 | ~1.007 |
| 1e16 | ~1.42 |
| 1e17 | ~1.64 |
Root cause: the transformed-rejection sampler random_poisson_ptrs
(numpy/random/src/distributions/distributions.c,
around L600) evaluates the log pmf: -lam + k*loglam - loggam(k+1), in
float64.
For example: for promising candidates (that is, those that should be accepted at
L600) k ≈ lam at lam = 1e17 the above three terms are ~1e17, ~3.91e18,
and ~3.81e18. They nearly cancel to an exact result of about ~ -20. But a
float64 computation returns 0.0 here due to rounding error, thereby never accepting the promising candidate k.
Reproduce the code example:
import numpy as np
lam = 1e17
size = 2_000_000
rng = np.random.default_rng()
x = rng.poisson(lam, size=size)
variance = x.var()
ratio = variance / lam
print(f"sample mean = {x.mean()}")
print(f"variance / lambda = {ratio}Error message:
No exception. The sampler silently returns samples from a wrong distribution.Python and NumPy Versions:
NumPy: 2.5.0
Python: 3.13.2
Runtime Environment:
[{'numpy_version': '2.5.0',
'python': '3.13.2 (main, Feb 4 2025, 14:51:09) [Clang 16.0.0 '
'(clang-1600.0.26.6)]',
'uname': uname_result(system='Darwin', node='infra04-wg313.visitor.eduroam.wireless.umn.edu', release='24.6.0', version='Darwin Kernel Version 24.6.0: Mon Jan 19 21:58:37 PST 2026; root:xnu-11417.140.69.708.3~1/RELEASE_ARM64_T8103', machine='arm64')},
{'simd_extensions': {'baseline': ['NEON', 'NEON_FP16', 'NEON_VFPV4', 'ASIMD'],
'found': ['ASIMDHP', 'ASIMDDP'],
'not_found': ['ASIMDFHM']}},
{'ignore_floating_point_errors_in_matmul': True}]
How does this issue affect you or how did you find it:
We use np.random.poisson as a backend in our model counting tool:
[pepin](https://github.com/meelgroup/pepin). For large model counts we received
wrong counts, which we traced back to this variance inflation in the Poisson
sampler at large lam.
Activity
WarrenWeckesser commentedon Jul 13, 2026
Thanks for reporting the issue and providing an analysis, @uddaloksarkar.
If we don't have a simple fix for the existing code, I suspect that for sufficiently large lambda we could use the normal distribution to approximate the Poisson distribution. For the values of lambda where the current Poisson algorithm starts to fail, the normal approximation should be pretty much indistinguishable from the true Poisson distribution.
adrian740 commentedon Jul 15, 2026
I agree with the analysis by @uddaloksarkar. I found the following resource: https://www.r-project.org/doc/reports/CLoader-dbinom-2002.pdf
At the end of page 3 and beginning of page 4, Loader documents the numerical instability of evaluating the direct Poisson log-PMF expression for large values of
lambdaandk. Equation 7 gives a numerically stable equivalent.Appendix A provides an implementation of
NumPy's PTRS acceptance test operates in log space, so the corresponding replacement for the right-hand side at L600 would be
Appendix A also provides implementations of both
stirlerrandbd0. Adapting these for NumPy should be relatively direct, with a separatek == 0case.Happy to open a PR implementing this if the approach sounds reasonable.
ngoldbaum commentedon Jul 16, 2026
Ping @rkern because this is an RNG internals issue.
rkern commentedon Jul 16, 2026
From my comment on the linked Scipy issue:
The paper @adrian740 found might just be such a one. It's unclear what the license of the code listings in the paper is as the linked website is long gone, but I don't think it's much of a problem as the text seems clear enough to implement from that.
rkern commentedon Jul 16, 2026
Ah, didn't scroll back far enough in the Wayback Machine. The code is permissively-licensed:
https://web.archive.org/web/20011227233957/http://cm.bell-labs.com/cm/ms/departments/sia/catherine/dbinom/dbinom.c