Direct Computation of 'ppois()' Poisson Distribution Probabilities

Direct computation and errors of ppois Poisson distribution probabilities.


ppoisD(q, lambda, all.from.0 = TRUE)
ppoisErr (lambda, ppFUN = ppoisD, iP = 1e-15,
          xM = qpois(iP, lambda=lambda, lower.tail=FALSE),
          verbose = FALSE)

numeric vector of non-negative integer values, “quantiles” at which to evaluate ppois(q, la) and ppFUN(q, la).


positive parameter of the Poisson distribution, lambda\(= \lambda = E[X] = Var[X]\) where \(X \sim Pois(\lambda)\).


logical indicating if q is positive integer, and the probabilities should computed for all quantile values of 0:q.


alternative ppois evaluation, by default the direct summation of dpois(k, lambda).


small number, iP << 1, used to construct the abscissa values x at which to evaluate and compare ppois() and ppFUN(), see xM:


(specified instead of iP: ) the maximal x-value to be used, i.e., the values used will be x <- 0:iM. The default, qpois(1-iP, lambda = lambda) is the upper tail iP-quantile of Poi(lambda).


logical indicating if extra information should be printed.


ppoisD() contains the poisson probabilities along q, i.e., is a numeric vector of length length(q).

re <- ppoisErr() returns the relative “error” of ppois(x0, lambda) where ppFUN(x0, lambda) is assumed to be the truth and x0 the “worst case”, i.e., the value (among x) with the largest such difference.

Additionally, attr(re, "x0") contains that value x0.

See Also


  • ppoisErr
  • ppoisD
(lams <- outer(c(1,2,5), 10^(0:3)))# 10^4 is already slow!
system.time(e1 <- sapply(lams, ppoisErr))
e1 / .Machine$double.eps

## Try another 'ppFUN' :---------------------------------
## this relies on the fact that it's *only* used on an 'x' of the form  0:M :
ppD0 <- function(x, lambda, all.from.0=TRUE)
            cumsum(dpois(if(all.from.0) 0:x else x, lambda=lambda))
## and test it:
p0 <- ppD0 (  1000, lambda=10)
p1 <- ppois(0:1000, lambda=10)
stopifnot(all.equal(p0,p1, tol=8*.Machine$double.eps))

system.time(p0.slow <- ppoisD(0:1000, lambda=10, all.from.0=FALSE))# not very slow, here
p0.1 <- ppoisD(1000, lambda=10)
if(requireNamespace("Rmpfr")) {
 ppoisMpfr <- function(x, lambda) cumsum(Rmpfr::dpois(x, lambda=lambda))
 p0.best <- ppoisMpfr(0:1000, lambda = Rmpfr::mpfr(10, precBits = 256))
 AllEq. <- Rmpfr::all.equal
 AllEq <- function(target, current, ...)
    AllEq.(target, current, ...,
           formatFUN = function(x, ...) Rmpfr::format(x, digits = 9))
 print(AllEq(p0.best, p0,      tol = 0)) # 2.06e-18
 print(AllEq(p0.best, p0.slow, tol = 0)) # the "worst" (4.44e-17)
 print(AllEq(p0.best, p0.1,    tol = 0)) # 1.08e-18

## Now (with 'all.from.0 = TRUE',  it is fast too):
p15    <- ppoisErr(2^13)
p15.0. <- ppoisErr(2^13, ppFUN = ppD0)
c(p15, p15.0.) / .Machine$double.eps # on Lnx 64b, see (-10  2.5), then (-2 -2)

## lapply(), so you see "x0" values :
str(e0. <- lapply(lams, ppoisErr, ppFUN = ppD0))

## The first version [called 'err.lambd0()' for years] used simple  cumsum(dpois(..))
## NOTE: It is *stil* much faster, as it relies on special  x == 0:M  relation
## Author: Martin Maechler, Date:  1 Mar 2004, 17:40
e0 <- sapply(lams, function(lamb) ppoisErr(lamb, ppFUN = ppD0))
all.equal(e1, e0) # typically TRUE,  though small "random" differences:
cbind(e1, e0) * 2^53 # on Lnx 64b, seeing integer values in {-24, .., 33}
# }
Documentation reproduced from package DPQ, version 0.3-3, License: GPL (>= 2)

Community examples

Looks like there are no examples yet.