R Under development (unstable) (2024-08-21 r87038 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > > 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) List of 2 $ maximum : num 1083.28151617 $ objective: num 2.48134846004e-14 > ## 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) [1] 2.481348460037225e-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)) [1] 1071 > ## 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)) [1] 0.7764318711312523 > (p2 <- sum(dpois(0:1001, lam))) [1] 0.7764318711312358 > 1 - p1/p2# 0, was -6.7e-6 [1] -2.131628207280301e-14 > stopifnot(abs(1 - p1/p2) < 1e-12) > > ## new pgamma(alphalimit = 1300): > lam <- 1274.487 > (p1 <- ppois(1301, lam)) [1] 0.775975173732999 > (p2 <- sum(dpois(0:1301, lam))) [1] 0.7759751737330141 > 1 - p1/p2# -2.22e-16, was -5.13e-6 [1] 1.942890293094024e-14 > 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) ) + ) [1] 0.775975173732999 > > > ##- 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()) [1] FALSE > > ## 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)) [1] 2.220446049250313e-16 attr(,"x") x0 xM 907 1148 > (pEd <- ppoisErr(900, ppFUN = ppD0)) [1] 2.220446049250313e-16 attr(,"x") x0 xM 904 1148 > options(digits = 7)# back to normal > all.equal(pE, pEd) # not quite: 'x0' differs slightly [1] "Attributes: < Component \"x\": Mean relative difference: 0.003307607 >" > > ## 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) + }) + } Lam = 11400; M = 11850 : Lam = 11401; M = 11850 : Lam = 11402.8; M = 11850 : Lam = 11408.1; M = 11850 : Lam = 11408.7; M = 11850 : Lam = 11415.5; M = 11850 : Lam = 11429.1; M = 11900 : Lam = 11431.7; M = 11900 : Lam = 11444.1; M = 11900 : Lam = 11528.5; M = 12000 : Lam = 12679.6; M = 13150 : > > > > ## 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")) [1] "D:/RCompile/CRANincoming/R-devel/lib/DPQ/safe" > 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) + } load(D:/RCompile/CRANincoming/R-devel/lib/DPQ/safe/tests_ppoisErr.rda) creates lams errL alphLim > > 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") #{lams > alphLim} : 0 > ## > 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----------------------------- > > > proc.time() user system elapsed 1.29 0.10 1.39