options(digits=16) Aeq <- function(x,y) all.equal(x,y, tol = 1e-14) stopifnot(exprs = { identical(ppois(1.5, 2), ppois(1, 2)) # by internal code definition Aeq(ppois(1,2), sum(dpois(0:1, 2))) }) ## From: Matthias Kohl ## To: r-help ## Subject: [R] Exactness of ppois ## Date: Thu, 15 Jan 2004 13:55:22 +0000 ##- by checking the precision of a convolution algorithm, we found the ##- following "inexactness": ##- We work with R Version 1.8.1 (2003-11-21) on Windows systems (NT, 2000, XP) ## Try the code: ## Kolmogorov distance between two methods to ## determine P(Poisson(lambda)<=x) Kolm.dist <- function(lam, eps) { x <- 0:qpois(eps, lambda = lam, lower.tail=FALSE) max(abs(ppois(x, lambda = lam) - cumsum(dpois(x, lambda = lam)))) } str(erg <- optimize(Kolm.dist, lower = 800, upper = 2000, maximum = TRUE, eps = 1e-15), digits=12) ## max at lambda = 1262.6, but value is simply 2.22e-16 ... ## i.e. very small! ==> no problem anymore Kolm.dist(lam = erg$max, eps = 1e-14) Kolm1.dist <- function(lam, eps) { x <- 0:qpois(eps, lambda = lam, lower.tail=FALSE) which.max(abs(ppois(x, lambda = lam) - cumsum(dpois(x, lambda = lam)))) } (x. <- Kolm1.dist(lam = erg$max, eps = 1e-15)) ## 1360, was 1422 # was '1001' ## or '1301' or .. : in any case == "alphlimit + 1" !! ##- So for lambda=977.8 and x=1001 we get a distance of about 5.2e-06. ##- (This inexactness seems to hold for all lambda values greater than ##- about 900.) ## MM: lam <- 977.8 (p1 <- ppois(1001, lam)) (p2 <- sum(dpois(0:1001, lam))) 1 - p1/p2# 0, was -6.7e-6 stopifnot(abs(1 - p1/p2) < 1e-12) ## new pgamma(alphalimit = 1300): lam <- 1274.487 (p1 <- ppois(1301, lam)) (p2 <- sum(dpois(0:1301, lam))) 1 - p1/p2# -2.22e-16, was -5.13e-6 stopifnot(abs(1 - p1/p2) < 1e-12) ## Also: lam <- 1274.487 x <- 1301 ## This is by R's definition of ppois() : stopifnot( pgamma(lam, x+1, lower=FALSE) == print( ppois(x,lam) ) ) ##- BUT, summing about 1000 terms of exactness around 1e-16, ##- we would expect an error of order 1e-13. ##- We suspect algorithm AS 239 (pgamma) to cause that flaw. ##- Do you think this could cause other problems apart from ##- that admittedly extreme example? ##- Thanks for your attention! ##- Matthias ##--- More generally: library(DPQ) (doExtras <- DPQ:::doExtras()) ## dyn.load("/u/maechler/R/MM/NUMERICS/dpq-functions/ppois-direct.so") ## "TODO": currently quite a few checks already happen in the examples in ## >> ../man/ppoisson.Rd << ppD0 <- function(x, lambda, all.from.0 = TRUE) cumsum(dpois(if(all.from.0) 0:x else x, lambda=lambda)) (pE <- ppoisErr(900)) (pEd <- ppoisErr(900, ppFUN = ppD0)) options(digits = 7)# back to normal all.equal(pE, pEd) # not quite: 'x0' differs slightly ## Large Lam where exp(-lambda) underflows to 0 (even in 'long double'): ## "the minimal" large lambda is okLD <- okLongDouble(999) notXorLD <- !okLD || !doExtras Lam <- 11400 ## but we choose a slightly larger, where one sees more set.seed(12) for(Lam in c(11400 + c(0, sort(round(rlnorm(10, 4, 2), 1))))) { M <- ceiling((Lam + 4*sqrt(Lam))/50)*50 tit <- sprintf("Lam = %g; M = %7.0f", Lam, M) cat(tit, ":\n") kM <- 0:M pp. <- ppoisD(M, lambda=Lam) ## with current DEBUG output vives ## ppoisD(*, lambda=11600): expl(-ldlam)=0= 0 ==> llam=9.35876, exp_arg=-11600 ## .. i=30, finally new f = expl(exp_arg = -11393.9) = 4.95747e-4949 > 0 ## ____________or________________ ## ppoisD(lambda=11444, expl(-ldlam)=0= 0 ==> llam=9.34522, exp_arg=-11444 ## .. i=6, finally new f = expl(exp_arg = -11394.5) = 2.69745e-4949 > 0 pp <- ppois(kM, lambda=Lam) pp0 <- cumsum(dp <- dpois(kM, lambda=Lam)) if(doExtras) print(system.time( ppSL <- ppoisD(kM, lambda=Lam, all.from.0=FALSE) )) ip <- dp > 0 plot(kM[ip], pp[ip], type="l", main = tit) plot(kM[ip], pp[ip], type="l", main = tit, log = "y")# LOG: here *no* warnings lines(kM[ip], pp0[ip], col=2) lines(kM[ip], pp.[ip], col=adjustcolor("blue", 1/2), lwd=4) abline(v = Lam, lty=2, col="gray") ## plot (kM[ip], pp.[ip] - pp[ip], xlab = "k", type="l", main = tit) lines(kM[ip], pp0[ip] - pp[ip], col = adjustcolor("red", 1/2), lwd = 2) if(doExtras) lines(kM[ip], ppSL[ip] - pp[ip], col = adjustcolor(3, 1/2), lwd=3) abline(v = Lam, lty=2, col="gray") stopifnot(exprs = if(doExtras) { !okLD || all.equal(pp, pp. , tol = 1e-12) notXorLD || all.equal(pp, ppSL, tol = 1e-12) notXorLD || all.equal(pp., ppSL, tol = 1e-14)# both ppoisD() "fast" or "slow" should be very close } else { !okLD || all.equal(pp, pp. , tol = 1e-12) }) } ## MM(2018-08): 'alphLim' below must have been a way to set 'alphlimit' in pgamma()'s C code. ## ---------- ## ## alphLim <- 1000 # old alphLim <- 100000 # new --- in R's C code for pgamma() since R 1.8.0 ## when using ppoisErr <- ppoisErr1 {much slower ==> save here: (sdir <- system.file("safe", package="DPQ")) sfil <- file.path(sdir, "tests_ppoisErr.rda") if(!doExtras && file.exists(sfil)) { cat(sprintf("load(%s) creates ", sfil), load(sfil), "\n") } else { ## 2^20 is much too large and the last few take much time! l2ex <- if(doExtras) seq(1, 15, by=1/8) else seq(1, 12, by=1/4) errL <- lams <- 2^l2ex ## Use for() loop instead of *apply() to see progress and time of each: cat(head <- paste( " lambda | errLam | user system elapsed [x 1000, i.e. ms]", "---------+------------+---------------------\n", sep="\n")) for(i in seq_along(lams)) { st <- system.time(errL[i] <- ppoisErr(lams[i])) cat(sprintf("%8.2f | %10.4g |", lams[i], errL[i]), paste(sprintf("%5.0f", 1000*as.vector(st)[1:3]), collapse=" "), "\n") } ; cat(head) save(lams,errL,alphLim, file = sfil) } plot(lams, abs(errL), log = "xy", xaxt="n", xlab = "lambda", main = paste("|rel.Error { ppois(x,lambda) } | (alphlimit=", formatC(alphLim)," in pgamma)")) abline(v= alphLim, lty=3) at.l <- c(outer(c(1,2,5),10^pretty(log10(lams)))) at.l <- at.l[10^par("usr")[1] < at.l & at.l < 10^par("usr")[2]] ## Improve axis(1, at=at.l, las=2) : axis(1, at=at.l, labels=NA) u.y <- par("usr")[3:4] yA <- function(e) 10^c(c(1+e, -e) %*% u.y) for(ll in at.l) { if(ll < 1e4) { # normal yt <- yA(.03); srt <- 0 ; adj <- 0.5 } else { # sloping yt <- yA(.03); srt <- -30 ; adj <- c(0.25, 0.5) } text(ll, yt, formatC(ll), srt = srt, adj = adj, xpd = NA) } lm1 <- lm(log10(abs(errL)) ~ log10(lams), subset = errL != 0 & lams < alphLim - 150) ## # 855.1, 975.5 already bad ^^^^^ abline(lm1, col = "blue") ## ----------------------- cat("#{lams > alphLim} : ", sum(lams > alphLim),"\n") ## if(sum(lams > alphLim) >= 2) { ## always FALSE nowadays lm2 <- lm(log10(abs(errL)) ~ log10(lams), subset = lams > alphLim) ## # alphlimit : ^^^^ abline(lm2, col = "purple") ## where do the lines cross: c1 <- coef(lm1) c2 <- coef(lm2) x. <- 10^print(log10x <- (c1[1]-c2[1])/(c2[2]-c1[2])) x. # 684'300 ## draw segment to cross point: points(x., 10^predict(lm1, new=data.frame(lams=x.)), type = 'h', col=2) text(x., 1e-16, paste("lambda=",formatC(x.)), srt=90, adj = c(0,0), col = 2) } # cannot draw these here-----------------------------