# NOT RUN {
(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}
# }
Run the code above in your browser using DataLab