Simple approximation of Poisson cumulative distribution in long tail?

I want to decide the capacity $C$ of a table so that it has residual odds less than $2^{-p}$ to overflow for given $p\in[40\dots 120]$, assuming the number of entries follows a Poisson law with a given expectancy $E\in[10^3\dots 10^{12}]$.

Ideally, I want the lowest integer C such that 1-CDF[PoissonDistribution[E],C] < 2^-p for given p and E; but I'm content with some C slightly higher than that. Mathematica is fine for manual computation, but I would like to compute C from p and E at compile time, which limits me to 64-bit integer arithmetic.

Update: In Mathematica (version 7) e = 1000; p = 40; c = Quantile[PoissonDistribution[e], 1 - 2^-p] is 1231 and seems about right (thanks @Procrastinator); however the result for both p = 50 and p = 60 is 1250, which is wrong on the unsafe side (and matters: my experiment repeats like $2^{25}$ times or more, and I want demonstrably less than $2^{-30}$ overall odds of failure). I want some crude but safe approximation using 64-bit integer arithmetic only, as available in C(++) at compile time.

• How about C = Quantile[PoissonDistribution[E],1-2^p]? – user10525 Sep 4 '12 at 9:50
• The leading term of the probability mass function of the Poisson dominates in the tail. – cardinal Sep 4 '12 at 10:25
• @Procrastinator: yes that works in Mathematica (except for sign of p, and precision issues, and names E and C that are reserved). BUT I need a simple approximation of that, possibly crude (but on the safe side) using 64-bit integer arityhmetic only! – fgrieu Sep 4 '12 at 10:37
• Re the update: Mathematica 8 returns 1262 for $p=50$ and 1290 for $p=60$. Re Normal approximation (@Proc): this cannot be expected to work well in the tails, which is crucial to the calculation. – whuber Sep 4 '12 at 11:41
• Perhaps you should ask on stackoverflow. I'm not familiar with the constraints you have. I don't know what stops you from using dynamic memory allocation, or whether you can use branching to decide the size of the array, or what the costs are of defining an array which is twice the size you need (and then not using all of it). If some function like $\mu + \sqrt{\log\log \mu} \log \mu \sqrt \mu + p \frac{\sqrt{\mu}} {\log \mu}$ (just as an example) gave you the exact answer, would you be able to implement an approximation under your constraints or not? It seems like a programming problem now. – Douglas Zare Sep 7 '12 at 0:26

A Poisson distribution with large mean is approximately normal, but you have to be careful that you want a tail bound and the normal approximation is proportionally less accurate near the tails.

One approach used in this MO question and with binomial distributions is to recognize that the tail decreases more rapidly than a geometric series, so you can write an explicit upper bound as a geometric series.

$$\begin{eqnarray}\sum_{k=D}^\infty \exp(-\mu)\frac{\mu^k}{k!} & \lt & \sum_{k=D}^\infty \exp(-\mu) \frac{\mu^D}{D!}\bigg(\frac \mu{D+1}\bigg)^{k-D} \\ & = & \exp(-\mu)\frac{\mu^D}{D!}\frac{1}{1-\frac{\mu}{D+1}} \\ & \lt & \exp(-\mu) \frac{\mu^D}{\sqrt{2\pi D}(D/e)^D} \frac{1}{1-\frac{\mu}{D+1}} \\ & = & \exp(D-\mu) \bigg(\frac{\mu}{D}\bigg)^D \frac{D+1}{\sqrt{2\pi D} (D+1-\mu)}\end{eqnarray}$$

Line 2 $\to$ line 3 was related to Stirling's formula. In practice I think you then want to solve $-p \log 2 = \log(\text{bound})$ numerically using binary search. Newton's method starting with an initial guess of $D = \mu + c \sqrt \mu.$ should also work.

For example, with $p=100$ and $\mu = 1000$, the numerical solution I get is 1384.89. A Poisson distribution with mean $1000$ takes the values from $0$ through $1384$ with probability $1-1/2^{100.06}.$ The values $0$ through $1383$ occur with probability $1-1/2^{99.59}.$

• +1. Another approach relates Poisson tail probabilities (on the right) to tail probabilities of Gamma distributions (on the left), which can be closely (over)estimated with a saddlepoint approximation. – whuber Sep 4 '12 at 12:20
• There's a long way from that to something restricted to 64-bit integer arithmetic (without exp, log, sqrt..) but I will work on it; thanks all! – fgrieu Sep 4 '12 at 14:15
• (+1) Up to the invocation of Stirling's approximation (which is irrelevant), this is exactly the bound I was (opaquely) referencing in my comment to the OP. (For example, see here.) – cardinal Sep 4 '12 at 20:16

You may see P. Harremoës: Sharp Bounds on Tail Probabilities for Poisson Random Variables https://helda.helsinki.fi/bitstream/handle/10138/229679/witmse_proc_17.pdf The main inequalities there are as follows. Let $$Y$$ be a Poisson random variable with parameter $$\lambda$$. Put $$G(x)= \sqrt{2\left(x\ln \frac{x}{\lambda} +\lambda-x\right)} \ \ {\rm sign} \left(x-\lambda\right).$$ Let $$\Phi$$ denote the cumulative distribution function for the standard normal law. Then, for all integer $$k\ge 0$$, $${\bf P}(Y which is equivalent to $$\Phi(G(k-1)) \le {\bf P}(Y Moreover, $$\Phi(G(k+(1/2))) \le {\bf P}(Y\le k).$$

• If you could write out the key equation (assuming there are only one or two) that would help in case the link goes dead at some time. – jbowman Nov 11 '18 at 3:00