# 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

$$\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}.$