Simple approximation of Poisson cumulative distribution in long tail?

Question

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.

Answer 1

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

Answer 2

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<k)\le \Phi(G(k)) \le {\bf P}(Y\le k),$$ which is equivalent to $$\Phi(G(k-1)) \le {\bf P}(Y<k)\le \Phi(G(k))$$ for all integer $k>0$. Moreover, $\Phi(G(k+(1/2))) \le {\bf P}(Y\le k)$ which implies that $$\Phi(G(k-1/2)) \le {\bf P}(Y<k)\le \Phi(G(k))$$ for all integer $k>0$.