Equation for Inverse Poisson CDF

Question

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?

Answer 1

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.