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. > #### qpois(), qbinom() and qnbinom() overflow much too early in the upper tail > #### ~~~~~ ~~~~~~ ~~~~~~~ > > if(!dev.interactive(orNone=TRUE)) pdf("qPoisBinom-ex.pdf") > .O.P. <- par(no.readonly=TRUE) > > ## Its source code, since {r80271 | 2021-05-08} in ~/R/D/r-devel/R/src/nmath/qDiscrete_search.h > ## contains > ' + /* Note : "same" code in qpois.c, qbinom.c, qnbinom.c -- + * FIXME: This is far from optimal [cancellation for p ~= 1, etc]: */ + if(!lower_tail || log_p) { + p = R_DT_qIv(p); /* need check again (cancellation!): */ + if (p == R_DT_0) return 0; + if (p == R_DT_1) return ML_POSINF; + } + /* temporary hack --- FIXME --- */ + if (p + 1.01*DBL_EPSILON >= 1.) return ML_POSINF; + ' [1] "\n /* Note : \"same\" code in qpois.c, qbinom.c, qnbinom.c --\n * FIXME: This is far from optimal [cancellation for p ~= 1, etc]: */\n if(!lower_tail || log_p) {\n\tp = R_DT_qIv(p); /* need check again (cancellation!): */\n\tif (p == R_DT_0) return 0;\n\tif (p == R_DT_1) return ML_POSINF;\n }\n /* temporary hack --- FIXME --- */\n if (p + 1.01*DBL_EPSILON >= 1.) return ML_POSINF;\n" > ## ==> same "bug" also in qbinom() and qnbinom() [! ?? ] > > e <- c(-1000, -500, -200, -100, seq(-70,-1/4, by=1/4)) > > lambda <- 10000 > qp <- qpois(2^e, lambda=lambda, lower.tail=FALSE) > > ## 'cbind_no_rownames' : > cbNoRN <- function(...) { + r <- cbind(...) + dimnames(r)[[1]] <- rep.int("", nrow(r)) + r + } > > cbNoRN(e, p=2^e, qpois=qp) e p qpois -1000.00 9.332636e-302 13934 -500.00 3.054936e-151 12728 -200.00 6.223015e-61 11687 -100.00 7.888609e-31 11170 -70.00 8.470329e-22 10967 -69.75 1.007298e-21 10965 -69.50 1.197885e-21 10963 -69.25 1.424534e-21 10961 -69.00 1.694066e-21 10960 -68.75 2.014595e-21 10958 -68.50 2.395771e-21 10956 -68.25 2.849068e-21 10954 -68.00 3.388132e-21 10952 -67.75 4.029190e-21 10950 -67.50 4.791542e-21 10948 -67.25 5.698136e-21 10946 -67.00 6.776264e-21 10945 -66.75 8.058381e-21 10943 -66.50 9.583084e-21 10941 -66.25 1.139627e-20 10939 -66.00 1.355253e-20 10937 -65.75 1.611676e-20 10935 -65.50 1.916617e-20 10933 -65.25 2.279254e-20 10931 -65.00 2.710505e-20 10929 -64.75 3.223352e-20 10927 -64.50 3.833234e-20 10925 -64.25 4.558509e-20 10923 -64.00 5.421011e-20 10921 -63.75 6.446705e-20 10920 -63.50 7.666467e-20 10918 -63.25 9.117017e-20 10916 -63.00 1.084202e-19 10914 -62.75 1.289341e-19 10912 -62.50 1.533293e-19 10910 -62.25 1.823403e-19 10908 -62.00 2.168404e-19 10906 -61.75 2.578682e-19 10904 -61.50 3.066587e-19 10902 -61.25 3.646807e-19 10900 -61.00 4.336809e-19 10898 -60.75 5.157364e-19 10896 -60.50 6.133174e-19 10894 -60.25 7.293614e-19 10892 -60.00 8.673617e-19 10890 -59.75 1.031473e-18 10888 -59.50 1.226635e-18 10886 -59.25 1.458723e-18 10884 -59.00 1.734723e-18 10882 -58.75 2.062946e-18 10880 -58.50 2.453269e-18 10878 -58.25 2.917446e-18 10876 -58.00 3.469447e-18 10874 -57.75 4.125891e-18 10872 -57.50 4.906539e-18 10870 -57.25 5.834891e-18 10868 -57.00 6.938894e-18 10866 -56.75 8.251782e-18 10863 -56.50 9.813078e-18 10861 -56.25 1.166978e-17 10859 -56.00 1.387779e-17 10857 -55.75 1.650356e-17 10855 -55.50 1.962616e-17 10853 -55.25 2.333956e-17 10851 -55.00 2.775558e-17 10849 -54.75 3.300713e-17 10847 -54.50 3.925231e-17 10845 -54.25 4.667913e-17 10843 -54.00 5.551115e-17 10840 -53.75 6.601426e-17 10838 -53.50 7.850462e-17 10836 -53.25 9.335826e-17 10834 -53.00 1.110223e-16 10832 -52.75 1.320285e-16 10830 -52.50 1.570092e-16 10828 -52.25 1.867165e-16 10826 -52.00 2.220446e-16 10823 -51.75 2.640570e-16 10821 -51.50 3.140185e-16 10819 -51.25 3.734330e-16 10817 -51.00 4.440892e-16 10815 -50.75 5.281140e-16 10812 -50.50 6.280370e-16 10810 -50.25 7.468660e-16 10808 -50.00 8.881784e-16 10806 -49.75 1.056228e-15 10804 -49.50 1.256074e-15 10802 -49.25 1.493732e-15 10799 -49.00 1.776357e-15 10797 -48.75 2.112456e-15 10795 -48.50 2.512148e-15 10793 -48.25 2.987464e-15 10790 -48.00 3.552714e-15 10788 -47.75 4.224912e-15 10786 -47.50 5.024296e-15 10784 -47.25 5.974928e-15 10781 -47.00 7.105427e-15 10779 -46.75 8.449825e-15 10777 -46.50 1.004859e-14 10775 -46.25 1.194986e-14 10772 -46.00 1.421085e-14 10770 -45.75 1.689965e-14 10768 -45.50 2.009718e-14 10765 -45.25 2.389971e-14 10763 -45.00 2.842171e-14 10761 -44.75 3.379930e-14 10758 -44.50 4.019437e-14 10756 -44.25 4.779943e-14 10754 -44.00 5.684342e-14 10751 -43.75 6.759860e-14 10749 -43.50 8.038873e-14 10747 -43.25 9.559885e-14 10744 -43.00 1.136868e-13 10742 -42.75 1.351972e-13 10740 -42.50 1.607775e-13 10737 -42.25 1.911977e-13 10735 -42.00 2.273737e-13 10732 -41.75 2.703944e-13 10730 -41.50 3.215549e-13 10728 -41.25 3.823954e-13 10725 -41.00 4.547474e-13 10723 -40.75 5.407888e-13 10720 -40.50 6.431099e-13 10718 -40.25 7.647908e-13 10715 -40.00 9.094947e-13 10713 -39.75 1.081578e-12 10710 -39.50 1.286220e-12 10708 -39.25 1.529582e-12 10705 -39.00 1.818989e-12 10703 -38.75 2.163155e-12 10700 -38.50 2.572439e-12 10698 -38.25 3.059163e-12 10695 -38.00 3.637979e-12 10693 -37.75 4.326310e-12 10690 -37.50 5.144879e-12 10688 -37.25 6.118327e-12 10685 -37.00 7.275958e-12 10683 -36.75 8.652621e-12 10680 -36.50 1.028976e-11 10677 -36.25 1.223665e-11 10675 -36.00 1.455192e-11 10672 -35.75 1.730524e-11 10670 -35.50 2.057952e-11 10667 -35.25 2.447331e-11 10664 -35.00 2.910383e-11 10662 -34.75 3.461048e-11 10659 -34.50 4.115903e-11 10656 -34.25 4.894661e-11 10654 -34.00 5.820766e-11 10651 -33.75 6.922096e-11 10648 -33.50 8.231806e-11 10646 -33.25 9.789323e-11 10643 -33.00 1.164153e-10 10640 -32.75 1.384419e-10 10638 -32.50 1.646361e-10 10635 -32.25 1.957865e-10 10632 -32.00 2.328306e-10 10629 -31.75 2.768839e-10 10627 -31.50 3.292723e-10 10624 -31.25 3.915729e-10 10621 -31.00 4.656613e-10 10618 -30.75 5.537677e-10 10615 -30.50 6.585445e-10 10612 -30.25 7.831458e-10 10610 -30.00 9.313226e-10 10607 -29.75 1.107535e-09 10604 -29.50 1.317089e-09 10601 -29.25 1.566292e-09 10598 -29.00 1.862645e-09 10595 -28.75 2.215071e-09 10592 -28.50 2.634178e-09 10589 -28.25 3.132583e-09 10586 -28.00 3.725290e-09 10583 -27.75 4.430142e-09 10580 -27.50 5.268356e-09 10577 -27.25 6.265167e-09 10574 -27.00 7.450581e-09 10571 -26.75 8.860283e-09 10568 -26.50 1.053671e-08 10565 -26.25 1.253033e-08 10562 -26.00 1.490116e-08 10559 -25.75 1.772057e-08 10556 -25.50 2.107342e-08 10553 -25.25 2.506067e-08 10550 -25.00 2.980232e-08 10547 -24.75 3.544113e-08 10544 -24.50 4.214685e-08 10540 -24.25 5.012133e-08 10537 -24.00 5.960464e-08 10534 -23.75 7.088227e-08 10531 -23.50 8.429370e-08 10527 -23.25 1.002427e-07 10524 -23.00 1.192093e-07 10521 -22.75 1.417645e-07 10518 -22.50 1.685874e-07 10514 -22.25 2.004853e-07 10511 -22.00 2.384186e-07 10508 -21.75 2.835291e-07 10504 -21.50 3.371748e-07 10501 -21.25 4.009707e-07 10497 -21.00 4.768372e-07 10494 -20.75 5.670581e-07 10490 -20.50 6.743496e-07 10487 -20.25 8.019413e-07 10483 -20.00 9.536743e-07 10480 -19.75 1.134116e-06 10476 -19.50 1.348699e-06 10473 -19.25 1.603883e-06 10469 -19.00 1.907349e-06 10466 -18.75 2.268233e-06 10462 -18.50 2.697398e-06 10458 -18.25 3.207765e-06 10454 -18.00 3.814697e-06 10451 -17.75 4.536465e-06 10447 -17.50 5.394797e-06 10443 -17.25 6.415531e-06 10439 -17.00 7.629395e-06 10435 -16.75 9.072930e-06 10432 -16.50 1.078959e-05 10428 -16.25 1.283106e-05 10424 -16.00 1.525879e-05 10420 -15.75 1.814586e-05 10416 -15.50 2.157919e-05 10412 -15.25 2.566212e-05 10408 -15.00 3.051758e-05 10403 -14.75 3.629172e-05 10399 -14.50 4.315837e-05 10395 -14.25 5.132424e-05 10391 -14.00 6.103516e-05 10386 -13.75 7.258344e-05 10382 -13.50 8.631675e-05 10378 -13.25 1.026485e-04 10373 -13.00 1.220703e-04 10369 -12.75 1.451669e-04 10364 -12.50 1.726335e-04 10360 -12.25 2.052970e-04 10355 -12.00 2.441406e-04 10351 -11.75 2.903338e-04 10346 -11.50 3.452670e-04 10341 -11.25 4.105940e-04 10336 -11.00 4.882812e-04 10331 -10.75 5.806675e-04 10326 -10.50 6.905340e-04 10321 -10.25 8.211879e-04 10316 -10.00 9.765625e-04 10311 -9.75 1.161335e-03 10306 -9.50 1.381068e-03 10301 -9.25 1.642376e-03 10295 -9.00 1.953125e-03 10290 -8.75 2.322670e-03 10284 -8.50 2.762136e-03 10279 -8.25 3.284752e-03 10273 -8.00 3.906250e-03 10267 -7.75 4.645340e-03 10261 -7.50 5.524272e-03 10255 -7.25 6.569503e-03 10249 -7.00 7.812500e-03 10243 -6.75 9.290681e-03 10236 -6.50 1.104854e-02 10230 -6.25 1.313901e-02 10223 -6.00 1.562500e-02 10216 -5.75 1.858136e-02 10209 -5.50 2.209709e-02 10202 -5.25 2.627801e-02 10194 -5.00 3.125000e-02 10187 -4.75 3.716272e-02 10179 -4.50 4.419417e-02 10171 -4.25 5.255603e-02 10162 -4.00 6.250000e-02 10154 -3.75 7.432544e-02 10145 -3.50 8.838835e-02 10135 -3.25 1.051121e-01 10125 -3.00 1.250000e-01 10115 -2.75 1.486509e-01 10104 -2.50 1.767767e-01 10093 -2.25 2.102241e-01 10081 -2.00 2.500000e-01 10067 -1.75 2.973018e-01 10053 -1.50 3.535534e-01 10037 -1.25 4.204482e-01 10020 -1.00 5.000000e-01 10000 -0.75 5.946036e-01 9976 -0.50 7.071068e-01 9945 -0.25 8.408964e-01 9900 > ## e p qpois > ## -1000.00 9.332636e-302 Inf > ## ..... ............ ... > ## ..... ............ ... > ## -52.25 1.867165e-16 Inf > ## -52.00 2.220446e-16 Inf > ## -51.75 2.640570e-16 Inf > ## -51.50 3.140185e-16 10770 > ## -51.25 3.734330e-16 10769 > ## -51.00 4.440892e-16 10769 > ## -50.75 5.281140e-16 10769 > ## -50.50 6.280370e-16 10769 > ## -50.25 7.468660e-16 10769 > ## -50.00 8.881784e-16 10769 > > plot(qp ~ e, type = "b", subset = -(1:5), + main = paste0("qpois(2^e, lambda=",lambda, + ") - early overflow to +Inf")) > > ## qbinom() "same" problem -- does not overflow to +Inf but to > ## 'size' ( = n in C ) = 100 here : as indeed, the source > ## ~/R/D/r-devel/R/src/nmath/qbinom.c has an *ADDITIONAL* hack here : > ' + q = 1 - pr; + if(q == 0.) return n; /* covers the full range of the distribution */ + ' [1] "\n q = 1 - pr;\n if(q == 0.) return n; /* covers the full range of the distribution */\n" > qBin <- qbinom(2^e, size = 100, prob = 0.4, lower.tail=FALSE) > cbNoRN(e, p=2^e, qBin)[c(1, 75:85),] e p qBin -1000.00 9.332636e-302 100 -52.50 1.570092e-16 80 -52.25 1.867165e-16 80 -52.00 2.220446e-16 80 -51.75 2.640570e-16 80 -51.50 3.140185e-16 80 -51.25 3.734330e-16 79 -51.00 4.440892e-16 79 -50.75 5.281140e-16 79 -50.50 6.280370e-16 79 -50.25 7.468660e-16 79 -50.00 8.881784e-16 79 > > plot(qBin ~ e, type = "b", subset = -(1:5), + main = paste0("qbinom(2^e, size = 100, prob = 0.4, lower.tail=FALSE", + ") - early overflow to 'size'")); abline(h=100, lty=3) > > ## qnbinom() "same" problem -- > ## ~/R/D/r-devel/R/src/nmath/qnbinom.c > qNB <- qnbinom(2^e, size = 100, prob = 0.4, lower.tail=FALSE) > cbNoRN(e, p=2^e, qNB)[c(1, 70:82),] e p qNB -1000.00 9.332636e-302 1949 -53.75 6.601426e-17 357 -53.50 7.850462e-17 357 -53.25 9.335826e-17 356 -53.00 1.110223e-16 355 -52.75 1.320285e-16 355 -52.50 1.570092e-16 354 -52.25 1.867165e-16 353 -52.00 2.220446e-16 353 -51.75 2.640570e-16 352 -51.50 3.140185e-16 351 -51.25 3.734330e-16 351 -51.00 4.440892e-16 350 -50.75 5.281140e-16 349 > ## e p qNB > ## -1000.00 9.332636e-302 0 > ## -53.75 6.601426e-17 0 > ## -53.50 7.850462e-17 0 <<<< !! even more wrong > ## -53.25 9.335826e-17 Inf > ## -53.00 1.110223e-16 Inf > ## -52.75 1.320285e-16 Inf > ## -52.50 1.570092e-16 Inf > ## -52.25 1.867165e-16 Inf > ## -52.00 2.220446e-16 Inf > ## -51.75 2.640570e-16 Inf << > ## -51.50 3.140185e-16 337 > ## -51.25 3.734330e-16 337 > ## -51.00 4.440892e-16 337 > ## -50.75 5.281140e-16 337 > > ## to make the jump to Inf, then 0, more visible, replace Inf by HUGE : > qN. <- qNB > qN.[qNB == Inf] <- 1e300 > plot(qN. ~ e, type = "l", subset = -(1:5), ylim = range(qNB, finite=TRUE), + main = paste0("qnbinom(2^e, size = 100, prob = 0.4, lower.tail=FALSE)", + ") - early \"overflow\" to **WRONG**")); abline(h=0, lty=3) > > require(DPQ) Loading required package: DPQ > ## Embarrassing forgotten FIXME (for boring boundary cases only): > M <- 2^31; pr <- 1e-9 > stopifnot(exprs = { + qbinomR (0:1, size=M, prob=pr) == c(0, M) # was 0 1 + qnbinomR(0:1, size=M, prob=pr) == c(0, Inf) # " " + qpoisR (0:1, M) == c(0, Inf) + qbinomR (c(-Inf,0), size=M, prob=pr, log.p=TRUE) == c(0, M) + qnbinomR(c(-Inf,0), size=M, prob=pr, log.p=TRUE) == c(0, Inf) + qpoisR (c(-Inf,0), M, log.p=TRUE) == c(0, Inf) + }) > > > > proc.time() user system elapsed 0.53 0.12 0.62