# Equation for Inverse Poisson CDF

I am attempting to calculate quantile probabilities. I.e., the value above which there is only a 1% chance occurrence for an arrival process.

The R code is pretty straight forward with say a lambda = 2. qpois(.99, 2, lower.tail = TRUE, log.p = FALSE)

But what is the inverse Poisson CDF equation?

• The help page for qpois explains the algorithm. It necessarily involves a numerical search.
– whuber
Aug 8, 2017 at 14:30
• @whuber That's interesting; I'd have assumed the quickest way would would be to use the relationship with the chi-square (and hence call an inverse incomplete gamma routine). Not having worked it through I'm probably missing something obvious though. Aug 9, 2017 at 3:16
• @glen The Gamma shape parameter is varied in the search.
– whuber
Aug 9, 2017 at 3:47
• @glen_b Your suggestion works pretty well. Here's a quick-and-dirty implementation. Half the code is just to find an upper bound on the solution. qpois0 <- Vectorize(function(q, lambda, ...) { f <- function(n) pgamma(lambda, n, lower.tail = FALSE, log.p = TRUE) - log(q); n.max <- pmax(1, qnorm(q, lambda, sqrt(lambda))); while(ppois(n.max, lambda) <= q) n.max <- 2 * n.max; floor(uniroot(f, c(0, n.max+1), ...)\$root) })
– whuber
Jul 13, 2022 at 16:04

Quantile function: If you already have the CDF $$F$$ available, you can write the quantile function for the Poisson distributon as:

\begin{align} Q(p) &\equiv \inf \Big\{ x = 0,1,2,... \Big| \ p \leqslant F(x) \Big\} \\[12pt] &= \sum_{x=0}^\infty \mathbb{I}(F(x) < p) \\[6pt] &= \sum_{x=0}^{x_U-1} \mathbb{I}(F(x) < p) \\[6pt] &= \sum_{x=0}^{x_U-1} \mathbb{I}(\log F(x) < \log p), \\[6pt] \end{align}

where $$x_U$$ is any non-negative integer that satisfies $$F(x_U) \geqslant p$$.$$^\dagger$$ I will call this formula the long-sum form. You can also write it as:

\begin{align} Q(p) &= x_L + \sum_{x=x_L}^{x_U-1} \mathbb{I}(\log F(x) < \log p), \ \ \ \\[6pt] \end{align}

where $$0 \leqslant x_L < x_U$$ is any non-negative integer that satisfies $$F(x_U) < p$$. I will call this latter formula the short-sum form. I have expressed both formulae in their final form using log-probabilities (instead of probabilities) because it is more accurate to use log-probabilities for computation.

Methods of computation: One simple way to construct an algorithm to compute the Poisson quantile function is to use the long-sum formula, with a first step that identifies an appropriate upper bound $$x_U$$ from some preliminary step. Alternatively, you can use the short-sum formula and search only over values between the lower and upper bound. Alternatively, you can use a "search method over the values between these bounds, where you directly look for the infimum value in the initial definition.

According to the documentation for the qpois function it uses the Cornish–Fisher Expansion followed by a search for the appropriate quantile value. The function calls the underlying C function C_qpois. The description suggests a search near the approximating value using the infimum definition.

Computation from the long-sum form: It is possible to compute the quantiles using the long-sum form by having a preliminary step that computes an upper bound $$M$$. There are a number of statistical papers that explore probability inequalities for the Poisson distribution (see e.g., Hoeffding 1963, Anderson and Sanders 1967, Short 2013). However, for simplicity, I will show how you can form a consevative bound using the Chebychev inequality. This inquality gives:

\begin{align} F(\lfloor \lambda + k\sqrt{\lambda} \rfloor) &\geqslant \mathbb{P}(\lambda - k\sqrt{\lambda} < X < \lambda + k\sqrt{\lambda}) > 1-\frac{1}{k^2}. \\[6pt] \end{align}

Setting $$p=1-1/k^2$$ gives $$k = 1/\sqrt{1-p}$$ so that:

$$x_U = \Bigg\lfloor \lambda + \sqrt{\frac{\lambda}{1-p}} \Bigg\rfloor.$$

We can implement this method in the function below. We treat the input probability $$p=1$$ as a special case giving the quantile $$q=\infty$$ as an output. We also allow inputs that are NA values and we replicate the output of the base function in this case. Remaining quantiles are computed from the long-sum formula above.

qpois.alt <- function(p, lambda, lower.tail = TRUE, log.p = FALSE) {

#Set log-probabilities (lower tail)
n <- length(p)
if (log.p) { LOGP <- p } else { LOGP <- log(p) }
if (!lower.tail) { LOGP <- VGAM::log1mexp(-LOGP) }

#Set output and deal with special cases (outputs NA and Inf)
QUANTILES <- as.numeric(rep(NA, n))
NNA   <- !is.na(LOGP)
NLOGP <- LOGP[NNA]
if (length(NLOGP) == 0) { return(QUANTILES) }
QUANTILES[NNA] <- rep(Inf, length(NNA))
if (min(NLOGP) >= 0)    { return(QUANTILES) }

#Set log-CDF vector
LP.MAX <- max(NLOGP[NLOGP < 0])
UPPER  <- floor(lambda + sqrt(lambda*exp(-VGAM::log1mexp(-LP.MAX))))
LOGCDF <- ppois(0:UPPER, lambda, lower.tail = TRUE, log.p = TRUE)

#Compute output
for (i in 1:n) {
if (NNA[i]) {
if (LOGP[i] < 0) {
QUANTILES[i] <- sum(LOGCDF < LOGP[i]) } } }

#Return output
QUANTILES }


Using the benchmark function it is possible to compare the run-speed of the alternative methods of programming the function. From the output below you will see that my version is a bit slower (but not much) than the version in the base package. (Note that the exact output from the benchmarking test depends on the speed of your computer and any background processes, etc., but the relative comparison should be reasonably stable over a large number of replications.)

#Benchmark the function against base version
library(rbenchmark)
PROBS <- 0:1000/1000
REPS  <- 10^5
BENCH <- benchmark('qpois'      = { qpois(PROBS, lambda = 2) },
'qpois.alt'  = { qpois.alt(PROBS, lambda = 2) },
replications = REPS)

#Show the benchmarking output
BENCH

test replications elapsed relative user.self sys.self user.child sys.child
1     qpois       100000   45.88    1.000     34.90    10.90         NA        NA
2 qpois.alt       100000   58.18    1.268     58.11     0.05         NA        NA


$$^\dagger$$ The bounds $$x_U$$ and $$x_L$$ implicitly depend on $$p$$, but we have written the notation without this dependence for brevity.