R Under development (unstable) (2023-11-02 r85465 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 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. > library(DPQ) > > ###---> Automatically find places where qgamma() is not so precise (PR#2214) : > ### For PR#2214, had '1e-8' below and found quite a bit > ## see /u/maechler/R/MM/NUMERICS/dpq-functions/beta-gamma-etc/qgamma-ex.R .. > > ## FIXME: Timing ! --- partly these matplot() partly get quite slow ~? > source(system.file(package="DPQ", "test-tools.R", mustWork=TRUE)) > ## => showProc.time(), ... list_() , loadList() , readRDS_() , save2RDS() > showProc.time() Time (user system elapsed): 0 0 0 > > (doExtras <- DPQ:::doExtras()) [1] FALSE > (sdir <- system.file("safe", package="DPQ")) ## save directory (to read from) [1] "D:/RCompile/CRANincoming/R-devel/lib/DPQ/safe" > > ### Nowadays finds cases in a special region for really small p and cutoff 1e-11 : > set.seed(47) > n <- if(doExtras) 100 else 32 > res <- cbind(p=1,df=1,rE=1)[-1,] > for(M in 1:(if(doExtras) 20 else 10)) + for(p in runif(n)) for(df in rlnorm(n)) { + r <- 1- pchisq(qchisq(p, df),df)/p + if(abs(r) > 1e-11) res <- rbind(res, c(p,df,r)) + } > > ### use df in U[0,1]: finds two cases with bound 1e-11 > for(p in runif(n)/2) for(df in runif(n)) { + qq <- qchisq(p, df) + if(qq > 0 && p > 0) { + r <- 1- pchisq(qq, df) / p + if(abs(r) > 1e-11) res <- rbind(res, c(p,df,r)) + } + } > > ### now df very close to 0 : ==> finds more cases > for(p in sort(c(runif(64)/2, exp(-(1+rlnorm(256)))))) + for(df in 2^-rlnorm(256, mean=2, sdlog=1.5)) { + qq <- qchisq(p, df) + if(qq > 0 && p > 0) { + r <- 1- pchisq(qq, df) / p + if(abs(r) > 1e-11) res <- rbind(res, c(p,df,r)) + } + } > showProc.time() Time (user system elapsed): 0.54 0.04 0.58 > > require(graphics) > if(!dev.interactive(orNone=TRUE)) pdf("qgamma-appr.pdf") > eaxis <- sfsmisc::eaxis > > showProc.time() Time (user system elapsed): 0.08 0 0.08 > ## if(nrow(res) > 0) { > cat("Found inaccurate examples where pchisq(qchisq(p, df),df) != p\n") Found inaccurate examples where pchisq(qchisq(p, df),df) != p > ## sort in p, then df: > res <- res[order(res[,"p"], res[,"df"]), ] > rE <- res[,"rE"] > if(nrow(res) > 20) { hist(rE, breaks = 30); rug(rE) } > plot(res[,1:2])##--> quite interesting : all along one curve > ## p <= 1/2 and df <= 1 (about) !! > res <- cbind(res, nDig = round(-log10(abs(rE)), 1)) > print(res, digits=12) p df rE nDig [1,] 0.000194375438651 0.02334079639198 -4.05718514340e-08 7.4 [2,] 0.000605300028912 0.02041606754775 -1.99857908001e-11 10.7 [3,] 0.001012316063255 0.01855615147677 -2.59145555106e-04 3.6 [4,] 0.001248285290785 0.01838201076117 -2.84196000067e-10 9.5 [5,] 0.001682388899865 0.01720736646288 5.53088974600e-04 3.3 [6,] 0.001746787400790 0.01731189518997 -5.86897839217e-08 7.2 [7,] 0.002664451237518 0.01599317398629 1.48421342013e-04 3.8 [8,] 0.002664451237518 0.01618024201222 -3.82806282229e-08 7.4 [9,] 0.003159421860255 0.01557612780310 -7.92117005632e-06 5.1 [10,] 0.003159421860255 0.01568183691729 -4.52237520765e-08 7.3 [11,] 0.004055462418244 0.01493858731306 4.15166391654e-06 5.4 [12,] 0.004400694140827 0.01459101672970 9.07907026434e-04 3.0 [13,] 0.004458811277768 0.01457506850867 9.03139988533e-05 4.0 [14,] 0.004481882165743 0.01468883074316 -3.23309491179e-07 6.5 [15,] 0.004939609905705 0.01440168350452 -2.81810098879e-06 5.6 [16,] 0.008824465120182 0.01276352706510 1.21107345756e-04 3.9 [17,] 0.009040265960535 0.01273711629661 1.38964402733e-05 4.9 [18,] 0.010839089634828 0.01242499920422 2.63413624246e-10 9.6 [19,] 0.011642124851282 0.01201471267173 1.44956234150e-04 3.8 [20,] 0.014753716559535 0.01155624353203 1.52962087441e-10 9.8 [21,] 0.015499213434879 0.01125420134457 -9.69695930770e-05 4.0 [22,] 0.015499213434879 0.01135920381800 -9.55739012376e-08 7.0 [23,] 0.018603016576955 0.01071716109330 1.63971046474e-03 2.8 [24,] 0.018603016576955 0.01073655493589 2.14388784340e-04 3.7 [25,] 0.022624242394389 0.01033379525113 -3.37865757594e-09 8.5 [26,] 0.022624242394389 0.01034206121729 -2.76332994265e-08 7.6 [27,] 0.023730217356634 0.01016252135853 -1.07732682708e-06 6.0 [28,] 0.032427027472295 0.00942923095016 5.11205522358e-11 10.3 [29,] 0.044753525441333 0.00839626444749 1.22224173549e-05 4.9 [30,] 0.081818424963746 0.00686007746204 8.92777740624e-10 9.0 [31,] 0.081818424963746 0.00689856335721 2.28502772259e-11 10.6 [32,] 0.082800309102258 0.00681234719059 4.17997558788e-09 8.4 [33,] 0.083507718914457 0.00680676700443 9.77167236016e-11 10.0 [34,] 0.090821658072474 0.00655269761981 -7.16033632386e-09 8.1 [35,] 0.102294760453517 0.00623563107239 3.69438657444e-09 8.4 [36,] 0.110869751789691 0.00603268830251 -3.44006823028e-10 9.5 [37,] 0.123950804624116 0.00571305309327 2.84683721041e-10 9.5 [38,] 0.127405857731893 0.00562369059572 6.60541454867e-09 8.2 [39,] 0.135229634154169 0.00540073357520 -2.34762594200e-05 4.6 [40,] 0.137732279982451 0.00533092076413 2.99285844990e-04 3.5 [41,] 0.138112917548194 0.00535138710974 -2.05335777981e-06 5.7 [42,] 0.141100635980184 0.00527305771429 4.31593832968e-05 4.4 [43,] 0.141100635980184 0.00537073537183 -3.00640179418e-10 9.5 [44,] 0.142905299416015 0.00523680041306 3.48180824883e-04 3.5 [45,] 0.145624557210331 0.00526923971034 -1.94501770245e-09 8.7 [46,] 0.154606872884529 0.00506806894407 -4.59924667240e-07 6.3 [47,] 0.154606872884529 0.00507366168703 2.72301046933e-07 6.6 [48,] 0.163535630067488 0.00497650928578 3.39664962823e-11 10.5 [49,] 0.169741036539408 0.00484181845356 5.31400978776e-09 8.3 [50,] 0.177327576288650 0.00465956102839 5.53404362603e-05 4.3 [51,] 0.178169157856761 0.00471949961255 4.79807527043e-10 9.3 [52,] 0.190094017358772 0.00450373552308 -1.29698447116e-06 5.9 [53,] 0.190147641510530 0.00453468705710 5.66235636157e-09 8.2 [54,] 0.200112534472267 0.00442273120514 7.20473680715e-11 10.1 [55,] 0.201518808589718 0.00439936964342 1.58748569845e-11 10.8 [56,] 0.201518808589718 0.00439976887947 -9.97182336704e-11 10.0 [57,] 0.210803673024037 0.00427351441034 -1.70232938856e-10 9.8 [58,] 0.213058614771766 0.00426179831847 1.10152997834e-11 11.0 [59,] 0.214780951412088 0.00419869272965 9.79194836326e-09 8.0 [60,] 0.232805106603566 0.00395399315002 -9.17581020055e-08 7.0 [61,] 0.249102914025652 0.00380019404026 -1.15818465929e-10 9.9 [62,] 0.249102914025652 0.00382493512126 -1.39670497390e-11 10.9 [63,] 0.252076511947811 0.00374903834738 -8.83337205604e-08 7.1 [64,] 0.253082914021191 0.00375259362798 3.65436092498e-09 8.4 [65,] 0.253922058700076 0.00371237348323 3.28994798726e-06 5.5 [66,] 0.254289278570932 0.00374343873151 -1.05664899053e-09 9.0 [67,] 0.260017499519858 0.00366179605930 2.34859742765e-07 6.6 [68,] 0.270323906831467 0.00351999192121 -1.56164756277e-04 3.8 [69,] 0.271699356057456 0.00355068132680 5.13092990317e-09 8.3 [70,] 0.275516196070002 0.00346804047756 -4.35171547588e-04 3.4 [71,] 0.280722231049885 0.00348224101220 5.48759926389e-10 9.3 [72,] 0.284601233201101 0.00344936339590 1.57145851887e-10 9.8 [73,] 0.290188543054775 0.00336613521112 -5.64443074502e-08 7.2 [74,] 0.290579022038283 0.00334423496113 1.02667567781e-07 7.0 [75,] 0.290579022038283 0.00336764858994 2.26061565023e-08 7.6 [76,] 0.291850198713803 0.00333552811650 -1.27338760580e-06 5.9 [77,] 0.296521136452775 0.00330308865102 2.25309977453e-07 6.6 [78,] 0.298034174946132 0.00330462333485 8.42470393447e-09 8.1 [79,] 0.300556783277253 0.00323922530004 4.66003314391e-05 4.3 [80,] 0.303182283998467 0.00328704590597 -1.46205270113e-11 10.8 [81,] 0.322319846303892 0.00306134512927 -1.15130830540e-05 4.9 [82,] 0.322319846303892 0.00310689001755 8.57751647487e-11 10.1 [83,] 0.325071272052651 0.00302343293053 -2.47088704493e-04 3.6 [84,] 0.325071272052651 0.00304146419577 3.18761056051e-06 5.5 [85,] 0.331888412404218 0.00300837121343 -4.96098895297e-09 8.3 [86,] 0.362278153188527 0.00278204202032 4.53939330569e-10 9.3 [87,] 0.385389476781711 0.00260981704384 7.37274796769e-10 9.1 [88,] 0.425333956955001 0.00232995789362 1.82823025607e-08 7.7 [89,] 0.439503709203564 0.00222452690840 -4.53585193982e-06 5.3 [90,] 0.439503709203564 0.00224964327069 -3.02331937263e-10 9.5 [91,] 0.450804624124430 0.00216770324934 -4.59455036239e-08 7.3 > > if(requireNamespace("scatterplot3d")) { + scatterplot3d::scatterplot3d(res[,1:3], type ='h') ## quite interesting: + ## the inaccurate (p,df) points are on nice monotone curve !!! + ## this is *less* revealing + scatterplot3d::scatterplot3d(res[,c("p","df","nDig")], type ='h') + } Loading required namespace: scatterplot3d > rL <- res[abs(res[,'rE']) > 1e-9,] > rL <- rL[order(rL[,1],rL[,2]),] > rL p df rE nDig [1,] 0.0001943754 0.023340796 -4.057185e-08 7.4 [2,] 0.0010123161 0.018556151 -2.591456e-04 3.6 [3,] 0.0016823889 0.017207366 5.530890e-04 3.3 [4,] 0.0017467874 0.017311895 -5.868978e-08 7.2 [5,] 0.0026644512 0.015993174 1.484213e-04 3.8 [6,] 0.0026644512 0.016180242 -3.828063e-08 7.4 [7,] 0.0031594219 0.015576128 -7.921170e-06 5.1 [8,] 0.0031594219 0.015681837 -4.522375e-08 7.3 [9,] 0.0040554624 0.014938587 4.151664e-06 5.4 [10,] 0.0044006941 0.014591017 9.079070e-04 3.0 [11,] 0.0044588113 0.014575069 9.031400e-05 4.0 [12,] 0.0044818822 0.014688831 -3.233095e-07 6.5 [13,] 0.0049396099 0.014401684 -2.818101e-06 5.6 [14,] 0.0088244651 0.012763527 1.211073e-04 3.9 [15,] 0.0090402660 0.012737116 1.389644e-05 4.9 [16,] 0.0116421249 0.012014713 1.449562e-04 3.8 [17,] 0.0154992134 0.011254201 -9.696959e-05 4.0 [18,] 0.0154992134 0.011359204 -9.557390e-08 7.0 [19,] 0.0186030166 0.010717161 1.639710e-03 2.8 [20,] 0.0186030166 0.010736555 2.143888e-04 3.7 [21,] 0.0226242424 0.010333795 -3.378658e-09 8.5 [22,] 0.0226242424 0.010342061 -2.763330e-08 7.6 [23,] 0.0237302174 0.010162521 -1.077327e-06 6.0 [24,] 0.0447535254 0.008396264 1.222242e-05 4.9 [25,] 0.0828003091 0.006812347 4.179976e-09 8.4 [26,] 0.0908216581 0.006552698 -7.160336e-09 8.1 [27,] 0.1022947605 0.006235631 3.694387e-09 8.4 [28,] 0.1274058577 0.005623691 6.605415e-09 8.2 [29,] 0.1352296342 0.005400734 -2.347626e-05 4.6 [30,] 0.1377322800 0.005330921 2.992858e-04 3.5 [31,] 0.1381129175 0.005351387 -2.053358e-06 5.7 [32,] 0.1411006360 0.005273058 4.315938e-05 4.4 [33,] 0.1429052994 0.005236800 3.481808e-04 3.5 [34,] 0.1456245572 0.005269240 -1.945018e-09 8.7 [35,] 0.1546068729 0.005068069 -4.599247e-07 6.3 [36,] 0.1546068729 0.005073662 2.723010e-07 6.6 [37,] 0.1697410365 0.004841818 5.314010e-09 8.3 [38,] 0.1773275763 0.004659561 5.534044e-05 4.3 [39,] 0.1900940174 0.004503736 -1.296984e-06 5.9 [40,] 0.1901476415 0.004534687 5.662356e-09 8.2 [41,] 0.2147809514 0.004198693 9.791948e-09 8.0 [42,] 0.2328051066 0.003953993 -9.175810e-08 7.0 [43,] 0.2520765119 0.003749038 -8.833372e-08 7.1 [44,] 0.2530829140 0.003752594 3.654361e-09 8.4 [45,] 0.2539220587 0.003712373 3.289948e-06 5.5 [46,] 0.2542892786 0.003743439 -1.056649e-09 9.0 [47,] 0.2600174995 0.003661796 2.348597e-07 6.6 [48,] 0.2703239068 0.003519992 -1.561648e-04 3.8 [49,] 0.2716993561 0.003550681 5.130930e-09 8.3 [50,] 0.2755161961 0.003468040 -4.351715e-04 3.4 [51,] 0.2901885431 0.003366135 -5.644431e-08 7.2 [52,] 0.2905790220 0.003344235 1.026676e-07 7.0 [53,] 0.2905790220 0.003367649 2.260616e-08 7.6 [54,] 0.2918501987 0.003335528 -1.273388e-06 5.9 [55,] 0.2965211365 0.003303089 2.253100e-07 6.6 [56,] 0.2980341749 0.003304623 8.424704e-09 8.1 [57,] 0.3005567833 0.003239225 4.660033e-05 4.3 [58,] 0.3223198463 0.003061345 -1.151308e-05 4.9 [59,] 0.3250712721 0.003023433 -2.470887e-04 3.6 [60,] 0.3250712721 0.003041464 3.187611e-06 5.5 [61,] 0.3318884124 0.003008371 -4.960989e-09 8.3 [62,] 0.4253339570 0.002329958 1.828230e-08 7.7 [63,] 0.4395037092 0.002224527 -4.535852e-06 5.3 [64,] 0.4508046241 0.002167703 -4.594550e-08 7.3 > plot(rL[,1:2], type = "b", main = "inaccurate pchisq/qchisq pairs") > > plot(rL[,1:2], type = "b", log = "x", ylim = range(0, rL[,"df"]), + xaxt = "n", + main = "inaccurate pchisq/qchisq pairs"); abline(h = 0, lty=2) > ## aha -- a perfect line !! > lines(res[,1:2], col = adjustcolor(1, 0.5)) > eaxis(1); axis(1, at = 1/2) > > d <- as.data.frame(res) > plot (df ~ log(p), data = d, type = "b", cex=1/4, col="gray") > points(df ~ log(p), data = as.data.frame(rL), col=2, cex = 1/2) > > summary(fm <- lm (df ~ log(p), data = d, weights = -log(abs(rE)))) Call: lm(formula = df ~ log(p), data = d, weights = -log(abs(rE))) Weighted Residuals: Min 1Q Median 3Q Max -6.924e-04 -1.443e-04 -2.096e-05 7.786e-05 1.079e-03 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 5.168e-06 1.149e-05 0.45 0.654 log(p) -2.725e-03 3.683e-06 -739.99 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.0002575 on 89 degrees of freedom Multiple R-squared: 0.9998, Adjusted R-squared: 0.9998 F-statistic: 5.476e+05 on 1 and 89 DF, p-value: < 2.2e-16 > ## R^2 = 0.9998 > > p0 <- 2^seq(-50,-1, by=1/8) > dN <- data.frame(p = p0, + df = predict(fm, newdata = data.frame(p = p0))) > rE <- with(dN, 1- pchisq(qchisq(p, df),df)/p) > dN <- cbind(dN, rE = rE, nDig = round(-log10(abs(rE)), 1)) > print(dN, digits=10) p df rE nDig 1 8.881784197e-16 0.094454797738 -6.103205421e-07 6.2 2 9.685654347e-16 0.094218673664 -2.417771916e-07 6.6 3 1.056228096e-15 0.093982549590 1.482102610e-07 6.8 4 1.151824906e-15 0.093746425517 5.596452862e-07 6.3 5 1.256073967e-15 0.093510301443 9.925313598e-07 6.0 6 1.369758374e-15 0.093274177369 -5.456116718e-07 6.3 7 1.493732098e-15 0.093038053295 -6.476817549e-08 7.2 8 1.628926404e-15 0.092801929221 4.375368099e-07 6.4 9 1.776356839e-15 0.092565805147 9.613067529e-07 6.0 10 1.937130869e-15 0.092329681073 -4.656781627e-07 6.3 11 2.112456192e-15 0.092093557000 1.060769860e-07 7.0 12 2.303649813e-15 0.091857432926 6.993075319e-07 6.2 13 2.512147934e-15 0.091621308852 -6.429917916e-07 6.2 14 2.739516748e-15 0.091385184778 -1.755251944e-09 8.8 15 2.987464197e-15 0.091149060704 6.609671205e-07 6.2 16 3.257852808e-15 0.090912936630 -5.966162631e-07 6.2 17 3.552713679e-15 0.090676812556 1.141328784e-07 6.9 18 3.874261739e-15 0.090440688483 8.463782823e-07 6.1 19 4.224912384e-15 0.090204564409 -3.264588258e-07 6.5 20 4.607299625e-15 0.089968440335 4.538341198e-07 6.3 21 5.024295868e-15 0.089732316261 -6.607810448e-07 6.2 22 5.479033495e-15 0.089496192187 1.675732128e-07 6.8 23 5.974928394e-15 0.089260068113 -8.888068794e-07 6.1 24 6.515705616e-15 0.089023944039 -1.237751479e-08 7.9 25 7.105427358e-15 0.088787819965 8.855724555e-07 6.1 26 7.748523477e-15 0.088551695892 -8.599116197e-08 7.1 27 8.449824769e-15 0.088315571818 8.600546403e-07 6.1 28 9.214599250e-15 0.088079447744 -5.324085839e-08 7.3 29 1.004859174e-14 0.087843323670 -9.348370564e-07 6.0 30 1.095806699e-14 0.087607199596 8.590026457e-08 7.1 31 1.194985679e-14 0.087371075522 -7.374084430e-07 6.1 32 1.303141123e-14 0.087134951448 3.314590533e-07 6.5 33 1.421085472e-14 0.086898827375 -4.335490789e-07 6.4 34 1.549704695e-14 0.086662703301 6.834623428e-07 6.2 35 1.689964954e-14 0.086426579227 -2.323215398e-08 7.6 36 1.842919850e-14 0.086190455153 -6.982055953e-07 6.2 37 2.009718347e-14 0.085954331079 4.935691396e-07 6.3 38 2.191613398e-14 0.085718207005 -1.230713460e-07 6.9 39 2.389971358e-14 0.085482082931 -7.079815811e-07 6.1 40 2.606282246e-14 0.085245958858 5.585871082e-07 6.3 41 2.842170943e-14 0.085009834784 3.202911969e-08 7.5 42 3.099409391e-14 0.084773710710 -4.627894383e-07 6.3 43 3.379929908e-14 0.084537586636 8.786038088e-07 6.1 44 3.685839700e-14 0.084301462562 4.421567634e-07 6.4 45 4.019436694e-14 0.084065338488 3.745831523e-08 7.4 46 4.383226796e-14 0.083829214414 -3.354886093e-07 6.5 47 4.779942715e-14 0.083593090340 -6.766810847e-07 6.2 48 5.212564492e-14 0.083356966267 7.928490540e-07 6.1 49 5.684341886e-14 0.083120842193 5.100597787e-07 6.3 50 6.198818782e-14 0.082884718119 2.590340996e-07 6.6 51 6.759859815e-14 0.082648594045 3.977492358e-08 7.4 52 7.371679400e-14 0.082412469971 -1.477148275e-07 6.8 53 8.038873388e-14 0.082176345897 -3.034322393e-07 6.5 54 8.766453592e-14 0.081940221823 -4.273743999e-07 6.4 55 9.559885430e-14 0.081704097750 -5.195383985e-07 6.3 56 1.042512898e-13 0.081467973676 -5.799213265e-07 6.2 57 1.136868377e-13 0.081231849602 -6.085202768e-07 6.2 58 1.239763756e-13 0.080995725528 -6.053323396e-07 6.2 59 1.351971963e-13 0.080759601454 -5.703546222e-07 6.2 60 1.474335880e-13 0.080523477380 -5.035842183e-07 6.3 61 1.607774678e-13 0.080287353306 -4.050182230e-07 6.4 62 1.753290718e-13 0.080051229233 -2.746537497e-07 6.6 63 1.911977086e-13 0.079815105159 -1.124878974e-07 6.9 64 2.085025797e-13 0.079578981085 8.148222741e-08 7.1 65 2.273736754e-13 0.079342857011 3.072595165e-07 6.5 66 2.479527513e-13 0.079106732937 5.648468603e-07 6.2 67 2.703943926e-13 0.078870608863 -8.276082293e-07 6.1 68 2.948671760e-13 0.078634484789 -5.012848909e-07 6.3 69 3.215549355e-13 0.078398360715 -1.431432306e-07 6.8 70 3.506581437e-13 0.078162236642 2.468196389e-07 6.6 71 3.823954172e-13 0.077926112568 6.686065978e-07 6.2 72 4.170051594e-13 0.077689988494 -5.340347400e-07 6.3 73 4.547473509e-13 0.077453864420 -4.348588734e-08 7.4 74 4.959055026e-13 0.077217740346 4.788953006e-07 6.3 75 5.407887852e-13 0.076981616272 -6.077622330e-07 6.2 76 5.897343520e-13 0.076745492198 -1.660406057e-08 7.8 77 6.431098711e-13 0.076509368125 6.063946640e-07 6.2 78 7.013162874e-13 0.076273244051 -3.642604891e-07 6.4 79 7.647908344e-13 0.076037119977 3.275302775e-07 6.5 80 8.340103188e-13 0.075800995903 -5.640569791e-07 6.2 81 9.094947018e-13 0.075564871829 1.965355126e-07 6.7 82 9.918110051e-13 0.075328747755 -6.159765531e-07 6.2 83 1.081577570e-12 0.075092623681 2.134273260e-07 6.7 84 1.179468704e-12 0.074856499608 -5.200022835e-07 6.3 85 1.286219742e-12 0.074620375534 3.782226260e-07 6.4 86 1.402632575e-12 0.074384251460 -2.761172810e-07 6.6 87 1.529581669e-12 0.074148127386 6.909382807e-07 6.2 88 1.668020638e-12 0.073912003312 1.156952842e-07 6.9 89 1.818989404e-12 0.073675879238 -4.174156929e-07 6.4 90 1.983622010e-12 0.073439755164 6.554522183e-07 6.2 91 2.163155141e-12 0.073203631090 2.014482300e-07 6.7 92 2.358937408e-12 0.072967507017 -2.104198267e-07 6.7 93 2.572439484e-12 0.072731382943 -5.801509226e-07 6.2 94 2.805265149e-12 0.072495258869 6.355293846e-07 6.2 95 3.059163338e-12 0.072259134795 3.449181628e-07 6.5 96 3.336041275e-12 0.072023010721 9.644778398e-08 7.0 97 3.637978807e-12 0.071786886647 -1.098807467e-07 7.0 98 3.967244020e-12 0.071550762573 -2.740664131e-07 6.6 99 4.326310282e-12 0.071314638500 -3.961082071e-07 6.4 100 4.717874816e-12 0.071078514426 -4.760051238e-07 6.3 101 5.144878969e-12 0.070842390352 -5.137561632e-07 6.3 102 5.610530299e-12 0.070606266278 -5.093603277e-07 6.3 103 6.118326675e-12 0.070370142204 -4.628166241e-07 6.3 104 6.672082550e-12 0.070134018130 -3.741240573e-07 6.4 105 7.275957614e-12 0.069897894056 -2.432816513e-07 6.6 106 7.934488041e-12 0.069661769983 -7.028841353e-08 7.2 107 8.652620563e-12 0.069425645909 1.448566248e-07 6.8 108 9.435749632e-12 0.069189521835 4.021544442e-07 6.4 109 1.028975794e-11 0.068953397761 7.016060164e-07 6.2 110 1.122106060e-11 0.068717273687 -4.176446862e-07 6.4 111 1.223665335e-11 0.068481149613 -2.873860860e-08 7.5 112 1.334416510e-11 0.068245025539 4.023232975e-07 6.4 113 1.455191523e-11 0.068008901466 -5.698258108e-07 6.2 114 1.586897608e-11 0.067772777392 -4.930793907e-08 7.3 115 1.730524113e-11 0.067536653318 5.133677967e-07 6.3 116 1.887149926e-11 0.067300529244 -3.116848801e-07 6.5 117 2.057951587e-11 0.067064405170 3.404481971e-07 6.5 118 2.244212120e-11 0.066828281096 -3.848110615e-07 6.4 119 2.447330670e-11 0.066592157022 3.567796256e-07 6.4 120 2.668833020e-11 0.066356032948 -2.686900891e-07 6.6 121 2.910383046e-11 0.066119908875 5.623584263e-07 6.2 122 3.173795216e-11 0.065883784801 3.667434378e-08 7.4 123 3.461048225e-11 0.065647660727 -4.365231998e-07 6.4 124 3.774299853e-11 0.065411536653 5.312784933e-07 6.3 125 4.115903175e-11 0.065175412579 1.578595065e-07 6.8 126 4.488424239e-11 0.064939288505 -1.630782611e-07 6.8 127 4.894661340e-11 0.064703164431 -4.315370075e-07 6.4 128 5.337666040e-11 0.064467040358 -6.475189360e-07 6.2 129 5.820766091e-11 0.064230916284 5.516172021e-07 6.3 130 6.347590433e-11 0.063994792210 4.354003288e-07 6.4 131 6.922096451e-11 0.063758668136 3.716548794e-07 6.4 132 7.548599706e-11 0.063522544062 3.603786357e-07 6.4 133 8.231806350e-11 0.063286419988 4.015693588e-07 6.4 134 8.976848478e-11 0.063050295914 4.952248247e-07 6.3 135 9.789322680e-11 0.062814171841 6.413427847e-07 6.2 136 1.067533208e-10 0.062578047767 -4.864749967e-07 6.3 137 1.164153218e-10 0.062341923693 -2.302655622e-07 6.6 138 1.269518087e-10 0.062105799619 7.839837712e-08 7.1 139 1.384419290e-10 0.061869675545 4.395145590e-07 6.4 140 1.509719941e-10 0.061633551471 -4.525762887e-07 6.3 141 1.646361270e-10 0.061397427397 1.860558396e-08 7.7 142 1.795369696e-10 0.061161303323 5.422316386e-07 6.3 143 1.957864536e-10 0.060925179250 -1.718041731e-07 6.8 144 2.135066416e-10 0.060689055176 4.618674431e-07 6.3 145 2.328306437e-10 0.060452931102 -1.317492586e-07 6.9 146 2.539036173e-10 0.060216807028 6.119535331e-07 6.2 147 2.768838580e-10 0.059980682954 1.387356798e-07 6.9 148 3.019439882e-10 0.059744558880 -2.716851295e-07 6.6 149 3.292722540e-10 0.059508434806 -6.193156095e-07 6.2 150 3.590739391e-10 0.059272310733 3.495619053e-07 6.5 151 3.915729072e-10 0.059036186659 1.222864675e-07 6.9 152 4.270132832e-10 0.058800062585 -4.221712024e-08 7.4 153 4.656612873e-10 0.058563938511 -1.439555917e-07 6.8 154 5.078072346e-10 0.058327814437 -1.829356824e-07 6.7 155 5.537677160e-10 0.058091690363 -1.591641428e-07 6.8 156 6.038879765e-10 0.057855566289 -7.264772095e-08 7.1 157 6.585445080e-10 0.057619442216 7.660682644e-08 7.1 158 7.181478783e-10 0.057383318142 2.885927388e-07 6.5 159 7.831458144e-10 0.057147194068 5.633032427e-07 6.2 160 8.540265665e-10 0.056911069994 -3.010137504e-07 6.5 161 9.313225746e-10 0.056674945920 1.043142719e-07 7.0 162 1.015614469e-09 0.056438821846 5.723448618e-07 6.2 163 1.107535432e-09 0.056202697772 -8.306449706e-08 7.1 164 1.207775953e-09 0.055966573698 5.155342467e-07 6.3 165 1.317089016e-09 0.055730449625 1.091102431e-09 9.0 166 1.436295757e-09 0.055494325551 -4.402707194e-07 6.4 167 1.566291629e-09 0.055258201477 3.567051224e-07 6.4 168 1.708053133e-09 0.055022077403 5.624482458e-08 7.2 169 1.862645149e-09 0.054785953329 -1.711695823e-07 6.8 170 2.031228938e-09 0.054549829255 -3.255506231e-07 6.5 171 2.215070864e-09 0.054313705181 -4.069108297e-07 6.4 172 2.415551906e-09 0.054077581108 -4.152627457e-07 6.4 173 2.634178032e-09 0.053841457034 -3.506189141e-07 6.5 174 2.872591513e-09 0.053605332960 -2.129918859e-07 6.7 175 3.132583258e-09 0.053369208886 -2.394214382e-09 8.6 176 3.416106266e-09 0.053133084812 2.811615332e-07 6.6 177 3.725290298e-09 0.052896960738 -4.754652252e-07 6.3 178 4.062457877e-09 0.052660836664 -4.082857119e-08 7.4 179 4.430141728e-09 0.052424712591 4.667263305e-07 6.3 180 4.831103812e-09 0.052188588517 -5.029512740e-08 7.3 181 5.268356064e-09 0.051952464443 -4.839869709e-07 6.3 182 5.745183026e-09 0.051716340369 2.526345069e-07 6.6 183 6.265166516e-09 0.051480216295 -1.969240593e-08 7.7 184 6.832212532e-09 0.051244092221 -2.087459210e-07 6.7 185 7.450580597e-09 0.051007968147 -3.145456327e-07 6.5 186 8.124915754e-09 0.050771844073 -3.371111472e-07 6.5 187 8.860283457e-09 0.050535720000 -2.764620672e-07 6.6 188 9.662207623e-09 0.050299595926 -1.326180012e-07 6.9 189 1.053671213e-08 0.050063471852 9.440144044e-08 7.0 190 1.149036605e-08 0.049827347778 4.045766356e-07 6.4 191 1.253033303e-08 0.049591223704 -2.420870762e-07 6.6 192 1.366442506e-08 0.049355099630 2.395534079e-07 6.6 193 1.490116119e-08 0.049118975556 -2.252219680e-07 6.6 194 1.624983151e-08 0.048882851483 4.277949085e-07 6.4 195 1.772056691e-08 0.048646727409 1.448079148e-07 6.8 196 1.932441525e-08 0.048410603335 -4.470065429e-08 7.3 197 2.107342426e-08 0.048174479261 -1.407587240e-07 6.9 198 2.298073210e-08 0.047938355187 -1.433942045e-07 6.8 199 2.506066606e-08 0.047702231113 -5.263500835e-08 7.3 200 2.732885013e-08 0.047466107039 1.314909438e-07 6.9 201 2.980232239e-08 0.047229982966 4.089557391e-07 6.4 202 3.249966302e-08 0.046993858892 -2.026429631e-07 6.7 203 3.544113383e-08 0.046757734818 2.666382147e-07 6.6 204 3.864883049e-08 0.046521610744 -1.427213290e-07 6.8 205 4.214684851e-08 0.046285486670 -4.483846485e-07 6.3 206 4.596146421e-08 0.046049362596 3.109955766e-07 6.5 207 5.012133212e-08 0.045813238522 2.073633717e-07 6.7 208 5.465770025e-08 0.045577114448 2.073183817e-07 6.7 209 5.960464478e-08 0.045340990375 3.108231487e-07 6.5 210 6.499932603e-08 0.045104866301 -4.225693508e-07 6.4 211 7.088226765e-08 0.044868742227 -1.068421185e-07 7.0 212 7.729766099e-08 0.044632618153 3.123191221e-07 6.5 213 8.429369702e-08 0.044396494079 -8.979923938e-08 7.0 214 9.192292842e-08 0.044160370005 -3.780712221e-07 6.4 215 1.002426642e-07 0.043924245931 3.616102581e-07 6.4 216 1.093154005e-07 0.043688121858 2.956315953e-07 6.5 217 1.192092896e-07 0.043451997784 3.433584006e-07 6.5 218 1.299986521e-07 0.043215873710 -3.936603832e-07 6.4 219 1.417645353e-07 0.042979749636 -1.134241190e-07 6.9 220 1.545953220e-07 0.042743625562 2.803691644e-07 6.6 221 1.685873940e-07 0.042507501488 -9.497757381e-08 7.0 222 1.838458568e-07 0.042271377414 -3.463649763e-07 6.5 223 2.004853285e-07 0.042035253341 3.982655141e-07 6.4 224 2.186308010e-07 0.041799129267 3.893573930e-07 6.4 225 2.384185791e-07 0.041563005193 -3.573736362e-07 6.4 226 2.599973041e-07 0.041326881119 -1.135260241e-07 6.9 227 2.835290706e-07 0.041090757045 2.539798885e-07 6.6 228 3.091906439e-07 0.040854632971 -1.007498536e-07 7.0 229 3.371747881e-07 0.040618508897 -3.214330786e-07 6.5 230 3.676917137e-07 0.040382384823 -4.081432381e-07 6.4 231 4.009706570e-07 0.040146260750 -3.609537462e-07 6.4 232 4.372616020e-07 0.039910136676 -1.799379763e-07 6.7 233 4.768371582e-07 0.039674012602 1.348307451e-07 6.9 234 5.199946082e-07 0.039437888528 -2.309643210e-07 6.6 235 5.670581412e-07 0.039201764454 3.563334906e-07 6.4 236 6.183812879e-07 0.038965640380 2.734302695e-07 6.6 237 6.743495762e-07 0.038729516306 3.344816787e-07 6.5 238 7.353834273e-07 0.038493392233 -2.537711121e-07 6.6 239 8.019413140e-07 0.038257268159 1.001817764e-07 7.0 240 8.745232040e-07 0.038021144085 -1.848126714e-07 6.7 241 9.536743164e-07 0.037785020011 -3.156868167e-07 6.5 242 1.039989216e-06 0.037548895937 -2.925440310e-07 6.5 243 1.134116282e-06 0.037312771863 -1.154876044e-07 6.9 244 1.236762576e-06 0.037076647789 2.153792595e-07 6.7 245 1.348699152e-06 0.036840523716 -5.632336997e-08 7.2 246 1.470766855e-06 0.036604399642 -1.638943072e-07 6.8 247 1.603882628e-06 0.036368275568 -1.074536522e-07 7.0 248 1.749046408e-06 0.036132151494 1.128786138e-07 6.9 249 1.907348633e-06 0.035896027420 -2.381848641e-07 6.6 250 2.079978433e-06 0.035659903346 3.148260224e-07 6.5 251 2.268232565e-06 0.035423779272 3.067288682e-07 6.5 252 2.473525152e-06 0.035187655198 -2.467805462e-07 6.6 253 2.697398305e-06 0.034951531125 9.809280876e-08 7.0 254 2.941533709e-06 0.034715407051 -9.217529584e-08 7.0 255 3.207765256e-06 0.034479282977 -9.840818915e-08 7.0 256 3.498092816e-06 0.034243158903 7.923719236e-08 7.1 257 3.814697266e-06 0.034007034829 -2.523280687e-07 6.6 258 4.159956866e-06 0.033770910755 2.978638364e-07 6.5 259 4.536465130e-06 0.033534786681 -3.332980627e-07 6.5 260 4.947050303e-06 0.033298662608 -8.305709964e-08 7.1 261 5.394796609e-06 0.033062538534 -3.110326443e-07 6.5 262 5.883067419e-06 0.032826414460 3.314533900e-07 6.5 263 6.415530512e-06 0.032590290386 -1.553353910e-07 6.8 264 6.996185632e-06 0.032354166312 2.279407841e-07 6.6 265 7.629394531e-06 0.032118042238 1.638951952e-07 6.8 266 8.319913732e-06 0.031881918164 3.135417902e-07 6.5 267 9.072930260e-06 0.031645794091 3.656239411e-08 7.4 268 9.894100606e-06 0.031409670017 -1.661171911e-08 7.8 269 1.078959322e-05 0.031173545943 1.537749448e-07 6.8 270 1.176613484e-05 0.030937421869 -7.675075087e-08 7.1 271 1.283106102e-05 0.030701297795 -7.365042554e-08 7.1 272 1.399237126e-05 0.030465173721 1.628068452e-07 6.8 273 1.525878906e-05 0.030229049647 2.398762633e-08 7.6 274 1.663982746e-05 0.029992925573 1.285369526e-07 6.9 275 1.814586052e-05 0.029756801500 -1.216179017e-07 6.9 276 1.978820121e-05 0.029520677426 -1.184344454e-07 6.9 277 2.157918644e-05 0.029284553352 1.377656036e-07 6.9 278 2.353226967e-05 0.029048429278 6.475007663e-08 7.2 279 2.566212205e-05 0.028812305204 2.546562354e-07 6.6 280 2.798474253e-05 0.028576181130 1.358057764e-07 6.9 281 3.051757812e-05 0.028340057056 -2.762893441e-07 6.6 282 3.327965493e-05 0.028103932983 1.553088368e-07 6.8 283 3.629172104e-05 0.027867808909 -2.519742930e-07 6.6 284 3.957640242e-05 0.027631684835 1.836182594e-07 6.7 285 4.315837288e-05 0.027395560761 -1.887623386e-07 6.7 286 4.706453935e-05 0.027159436687 -2.586992647e-07 6.6 287 5.132424410e-05 0.026923312613 -2.666509480e-08 7.6 288 5.596948506e-05 0.026687188539 -2.210355921e-08 7.7 289 6.103515625e-05 0.026451064466 -2.296375001e-07 6.6 290 6.655930986e-05 0.026214940392 -1.155410556e-07 6.9 291 7.258344208e-05 0.025978816318 -1.934344138e-07 6.7 292 7.915280485e-05 0.025742692244 5.977599926e-08 7.2 293 8.631674575e-05 0.025506568170 1.410128029e-07 6.9 294 9.412907870e-05 0.025270444096 6.554036214e-08 7.2 295 1.026484882e-04 0.025034320022 -1.513968051e-07 6.8 296 1.119389701e-04 0.024798195948 -7.984727768e-09 8.1 297 1.220703125e-04 0.024562071875 1.377294434e-08 7.9 298 1.331186197e-04 0.024325947801 -7.096806165e-08 7.1 299 1.451668842e-04 0.024089823727 2.236347701e-07 6.7 300 1.583056097e-04 0.023853699653 -3.406037208e-08 7.5 301 1.726334915e-04 0.023617575579 1.071465777e-07 7.0 302 1.882581574e-04 0.023381451505 1.915567783e-07 6.7 303 2.052969764e-04 0.023145327431 -2.153880070e-07 6.7 304 2.238779402e-04 0.022909203358 -1.943811714e-07 6.7 305 2.441406250e-04 0.022673079284 -1.853585185e-07 6.7 306 2.662372394e-04 0.022436955210 -1.734718695e-07 6.8 307 2.903337683e-04 0.022200831136 -1.439218822e-07 6.8 308 3.166112194e-04 0.021964707062 -8.196068868e-08 7.1 309 3.452669830e-04 0.021728582988 2.710549540e-08 7.6 310 3.765163148e-04 0.021492458914 1.979137895e-07 6.7 311 4.105939528e-04 0.021256334841 3.784054214e-08 7.4 312 4.477558805e-04 0.021020210767 -2.082511652e-08 7.7 313 4.882812500e-04 0.020784086693 3.635452528e-08 7.4 314 5.324744788e-04 0.020547962619 -1.675743373e-07 6.8 315 5.806675366e-04 0.020311838545 1.696075687e-07 6.8 316 6.332224388e-04 0.020075714471 -9.599202855e-08 7.0 317 6.905339660e-04 0.019839590397 -1.676104118e-07 6.8 318 7.530326296e-04 0.019603466323 -3.122729564e-08 7.5 319 8.211879055e-04 0.019367342250 -3.779823099e-08 7.4 320 8.955117609e-04 0.019131218176 -1.576477497e-07 6.8 321 9.765625000e-04 0.018895094102 -6.902103111e-09 8.2 322 1.064948958e-03 0.018658970028 7.899459398e-08 7.1 323 1.161335073e-03 0.018422845954 1.293463429e-07 6.9 324 1.266444878e-03 0.018186721880 -1.651599340e-07 6.8 325 1.381067932e-03 0.017950597806 -9.328011985e-08 7.0 326 1.506065259e-03 0.017714473733 3.008481086e-08 7.5 327 1.642375811e-03 0.017478349659 -8.903183302e-08 7.1 328 1.791023522e-03 0.017242225585 -8.893824854e-08 7.1 329 1.953125000e-03 0.017006101511 5.866337793e-08 7.2 330 2.129897915e-03 0.016769977437 7.502566035e-08 7.1 331 2.322670146e-03 0.016533853363 3.812256177e-09 8.4 332 2.532889755e-03 0.016297729289 -1.115639645e-07 7.0 333 2.762135864e-03 0.016061605216 6.307704159e-08 7.2 334 3.012130518e-03 0.015825481142 -1.675983974e-08 7.8 335 3.284751622e-03 0.015589357068 -1.223428825e-08 7.9 336 3.582047044e-03 0.015353232994 1.188536286e-07 6.9 337 3.906250000e-03 0.015117108920 -1.217006742e-07 6.9 338 4.259795831e-03 0.014880984846 -1.314145717e-07 6.9 339 4.645340293e-03 0.014644860772 -1.288061278e-07 6.9 340 5.065779510e-03 0.014408736698 -5.759581789e-08 7.2 341 5.524271728e-03 0.014172612625 -1.110605852e-07 7.0 342 6.024261037e-03 0.013936488551 2.552848366e-08 7.6 343 6.569503244e-03 0.013700364477 -7.038977312e-08 7.2 344 7.164094088e-03 0.013464240403 -8.014783304e-08 7.1 345 7.812500000e-03 0.013228116329 6.514375850e-08 7.2 346 8.519591661e-03 0.012991992255 -1.226700630e-08 7.9 347 9.290680586e-03 0.012755868181 3.829303208e-09 8.4 348 1.013155902e-02 0.012519744108 -1.720012133e-08 7.8 349 1.104854346e-02 0.012283620034 2.105846186e-08 7.7 350 1.204852207e-02 0.012047495960 1.170781927e-08 7.9 351 1.313900649e-02 0.011811371886 6.422402921e-08 7.2 352 1.432818818e-02 0.011575247812 9.486124963e-08 7.0 353 1.562500000e-02 0.011339123738 3.894983924e-08 7.4 354 1.703918332e-02 0.011102999664 3.208461541e-08 7.5 355 1.858136117e-02 0.010866875591 3.155170392e-08 7.5 356 2.026311804e-02 0.010630751517 1.297325769e-08 7.9 357 2.209708691e-02 0.010394627443 -3.001927462e-08 7.5 358 2.409704415e-02 0.010158503369 7.470776675e-08 7.1 359 2.627801298e-02 0.009922379295 2.881776118e-08 7.5 360 2.865637635e-02 0.009686255221 4.355808025e-08 7.4 361 3.125000000e-02 0.009450131147 3.486255895e-08 7.5 362 3.407836665e-02 0.009214007073 -4.540194531e-08 7.3 363 3.716272234e-02 0.008977883000 6.057274282e-08 7.2 364 4.052623608e-02 0.008741758926 -4.748040294e-08 7.3 365 4.419417382e-02 0.008505634852 -4.947996102e-08 7.3 366 4.819408829e-02 0.008269510778 5.303486095e-09 8.3 367 5.255602595e-02 0.008033386704 2.750289219e-09 8.6 368 5.731275270e-02 0.007797262630 7.649781364e-09 8.1 369 6.250000000e-02 0.007561138556 2.575243696e-08 7.6 370 6.815673329e-02 0.007325014483 2.588699921e-08 7.6 371 7.432544469e-02 0.007088890409 -3.873516863e-08 7.4 372 8.105247217e-02 0.006852766335 -2.868915505e-08 7.5 373 8.838834765e-02 0.006616642261 8.820975506e-09 8.1 374 9.638817659e-02 0.006380518187 -1.249981407e-08 7.9 375 1.051120519e-01 0.006144394113 -2.542299260e-08 7.6 376 1.146255054e-01 0.005908270039 -3.411815497e-08 7.5 377 1.250000000e-01 0.005672145966 1.097270552e-08 8.0 378 1.363134666e-01 0.005436021892 -5.883867171e-09 8.2 379 1.486508894e-01 0.005199897818 -1.778870562e-08 7.7 380 1.621049443e-01 0.004963773744 -2.477515548e-08 7.6 381 1.767766953e-01 0.004727649670 -2.316978920e-08 7.6 382 1.927763532e-01 0.004491525596 -1.391214588e-08 7.9 383 2.102241038e-01 0.004255401522 4.032062351e-09 8.4 384 2.292510108e-01 0.004019277448 1.322723042e-10 9.9 385 2.500000000e-01 0.003783153375 1.667484040e-09 8.8 386 2.726269332e-01 0.003547029301 -6.843556166e-09 8.2 387 2.973017788e-01 0.003310905227 -5.830291139e-09 8.2 388 3.242098887e-01 0.003074781153 -4.568762479e-09 8.3 389 3.535533906e-01 0.002838657079 -6.372726347e-09 8.2 390 3.855527064e-01 0.002602533005 4.293519318e-09 8.4 391 4.204482076e-01 0.002366408931 4.466780279e-09 8.4 392 4.585020216e-01 0.002130284858 -3.886752920e-09 8.4 393 5.000000000e-01 0.001894160784 1.495882751e-09 8.8 > > ## } ## only when we find inaccurate regions > showProc.time() Time (user system elapsed): 0.1 0.01 0.11 > > > ## Oops: another qgamma() / qchisq() problem: mostly NaN's == all solved now > curve(qgamma(x, 20), 1e-16, 1e-10, log='x') > curve(qgamma(x, 20), 1e-300, .99 , log='xy') # and add the critical region from above: > abline(v=c(1e-16, 1e-10), col="light blue") > curve(qgamma(x, 20), 1e-26, 1e-07, log='x') > ##-> now using log=TRUE in same region: > curve(qgamma(x, 20, log=TRUE), -38, -16)## no problem!! > curve(qgamma(exp(x), 20), add=TRUE, col="green3", n=2001) > ## had problem here, but no longer ! > > ##--> Further fix for qgamma: when 'x' is very small: use "log=TRUE of log(x)"! > > ## had bug (gave NaN), but no longer: > (q_12 <- qgamma(1e-12, 20)) [1] 2.330042 > all.equal(1e-12, pgamma(q_12, 20), tol=0)# show rel.err (Lnx 64-bit: 4.04e-16) [1] "Mean relative difference: 4.038968e-16" > stopifnot( + all.equal(1e-12, pgamma(q_12, 20), tolerance = 1e-14) + ) > > > ## --- Nice graphic : --- but amazingly *S..L..O..W* > > p.qgammaSml <- function(from= 1e-110, to = 1e-5, ylim = c(0.4, 1000), + n = 201, k.lab = 3, + a1 = c(10, seq(10.1,20, by=.2), 21:105), + a2 = seq(110,330, by=10), + a3 = seq(350,1600, by=50)) + { + ## Purpose: nice qgamma() lines ``for small x'' aka p + ## ---------------------------------------------------------------------- + ## Author: Martin Maechler, Date: 22 Mar 2004, 14:23 + x <- exp(seq(log(from), log(to), length = n)) + + op <- par(las=1, lab = c(10,10, 7), xaxs = "i", mex = 0.8) + on.exit(par(op)) + plot(x, qgamma(x, a1[1]), log="xy", ylim=ylim, type='l', xaxt = "n", + main = paste("qgamma(x, a) for very small x, a in [", + formatC(a1[1]),", ",formatC(max(a1,a2,a3)),"] - log-log", sep=''), + sub = R.version.string) + lab.x <- pretty(log10(c(from,to)), 20) + axis(1, at=10^lab.x, lab = paste("10^",formatC(lab.x),sep='')) + if(is.nan(qgamma(1e-12, 20))) + text(1e-60, 20, "all NaN", cex = 2) + if(!is.finite(qgamma(1e-140, 155))) + text(1e-240, 5, "all +Inf", cex = 2) + + lines.txt <- function(a.s, col = par("col")) { + col <- rep(col, length=length(a.s)) + for(i in seq(along=a.s)) { + qx <- qgamma(x, (a <- a.s[i])) + if(i %% k.lab == 0 && + any(ifi <- is.finite(qx) & qx >= ylim[1])) { + ik <- (i%%(2*k.lab))/k.lab # = 0 or 1 + j <- quantile(which(ifi), c(.02,(1:3)/4+ ik/10, .98)) + ## "segments" around the labels : + i0 <- 1 + for(jj in j) { + ii <- i0:(jj-1) + i2 <- jj + -1:1 + lines(x[ii], qx[ii], col=col[i]) + lines(x[i2], qx[i2], col=col[i], type = 'c') + i0 <- jj+1 + } + text(x[j], qx[j], formatC(a), col= "gray40", cex = 0.8) + } + else + lines(x, qx, col=col[i]) + + } + } + oo <- options(warn = -1) + lines.txt(a1[-1]) + lines.txt(a2, col= 2) + lines.txt(a3, col= rainbow(length(a3), .8, .8, + start = (max(a3)-min(a3))/(1+max(a3)))) + invisible(options(oo)) + } > > showProc.time() Time (user system elapsed): 0.01 0.02 0.03 > > p.qgammaSml() > p.qgammaSml(1e-300) > p.qgammaSml(1e-300,1e-50, a2= seq(100,360, by=4), a3=seq(350,1500, by=10)) > > showProc.time() Time (user system elapsed): 0.99 0.08 1.06 > > ## The "upper" problematic corner: > p.qgammaSml(1e-19, 1e-3, a2=NULL,a3=NULL, ylim=c(.1,20)) > p.qgammaSml(1e-19, 1e-3, a2=seq(1,12, by=.04), ylim=c(.1,20),a3=NULL,k.lab=10) > ## now shows the problem (quite well): > ## could it be in pgamma()'s inaccuracy, leading to qgamma() bias ? > aa <- c(seq(9, 22, by=0.25),seq(22.3,40,by=0.4)) > caa <- formatC(range(aa)) > sfsmisc::mult.fig(2) > curve(pgamma(x, sh=aa[1]), 0.5, 20, log = 'xy', ylim = c(1e-60, .2), + main = sprintf("pgamma(x, a) for a in [%s,%s]", caa[1],caa[2])) > for(sh in aa) curve(pgamma(x, sh), add = TRUE, col=2) > abline(h=c(1e-15), col="light blue", lty=2) > > curve(pgamma(x, sh=aa[1]), 0.5, 20, log = 'xy', ylim = c(1e-15, .8), + main = sprintf("pgamma(x, a) for a in [%s,%s]", caa[1],caa[2])) > for(sh in aa) curve(pgamma(x, sh), add = TRUE, col=2) > ## the "border curve" between "Pearson" and "Continued fraction (upper tail)" > ## in pgamma.c : > curve(pgamma(max(1,x), x), add = TRUE, col=4) > ## ==> pgamma() is perfect here {series expansion up to eps_C accuracy}! > > aa <- c(seq(9, 22, by=0.25),seq(22.3,40.4,by=0.4)) > p.qgammaSml(1e-24, 1e-5, a1=aa, a2=NULL,a3=NULL, ylim=c(.8,8)) > ## -------- save the above? > aa1 <- c(aa,seq(40.5,90, by=0.5)) > p.qgammaSml(1e-60, 1e-5, a1=aa1, a2=NULL,a3=NULL, ylim=c(.9, 16)) > aa2 <- c(aa1, seq(91,150, by= 1)) > p.qgammaSml(1e-90, 1e-5, a1=aa2, a2=NULL,a3=NULL, ylim=c(.9, 35)) > aa3 <- c(aa2, seq(150,250, by= 2), seq(253, 400, by=5)) > p.qgammaSml(1e-200, 1e-5, a1=aa3, a2=NULL,a3=NULL, ylim=c(.9, 100)) > p.qgammaSml(1e-200, 1e-5, a1=aa3, a2=NULL,a3=NULL, ylim=c(.9, 200),k.lab=9e9) > p.qgammaSml(1e-60, 1e-5, a1=aa3, a2=NULL,a3=NULL, ylim=c(.9, 200),k.lab=9e9) > > showProc.time() Time (user system elapsed): 3.11 0.09 3.2 > > ## lower a \> 10 > > curve(qgamma(x, 19), 1e-14, 1e-9, log='x') > curve(qgamma(x, 18), 1e-14, 1e-9, log='x') > curve(qgamma(x, 15), 1e-11, 5e-9, log='x') > curve(qgamma(x, 13), 5e-10, 1e-8, log='x') > curve(qgamma(x, 11), 1e-8, 5e-8, log='x') > curve(qgamma(x, 10.5), 4.2e-8, 6e-8, log='x') > curve(qgamma(x, 10.3), 6e-8, 7e-8, log='x') > curve(qgamma(x, 10.2), 7.1e-8, 7.6e-8, log='x') > curve(qgamma(x, 10.15),7.7e-8, 7.9e-8, log='x') > curve(qgamma(x, 10.14),7.88e-8,7.92e-8, log='x',n=10001) > > ## no more problems for smaller a!! here: > curve(qgamma(x, 10.13), 1e-10, 5e-4, log='x',n=20001) > curve(qgamma(x, 10.12), 1e-10, 5e-4, log='x',n=20001) > curve(qgamma(x, 10.1), 1e-10, 5e-4, log='x',n=20001) > > showProc.time() Time (user system elapsed): 0.62 0 0.63 > > ##--- the "+Inf" / premature "0" case: > curve(qgamma(x, 155, log=TRUE), -1500, 0, log='y', n=2001,col=2) > curve(qgamma(x, 1e3, log=TRUE), -1500, 0, log='y', n=2001,col=2) > ## now works, but slowly and with kink > curve(qgamma (x, 1e5, log=TRUE), -3e5, 0, log='y', n=2001,col=2,lwd=3) > curve(qgammaAppr(x, 1e5, log=TRUE), add = TRUE, n=2001, col="blue",lwd=.4) > ## --- curves are almost "identical" > ## ===> the kink *does* come from the initial approx... hmm > > ## still "identical" > curve(qgamma (x, 1e4, log=TRUE), -3e4, 0, log='y', n=2001,col=2) > curve(qgammaAppr(x, 1e4, log=TRUE), add = TRUE, n=2001, col="tomato3") > > ## now see some difference (approx. has kink at ~ -165) > curve(qgamma (x, 100, log=TRUE), -200, 0, log='y', n=2001,col=2) > curve(qgammaAppr(x, 100, log=TRUE), add = TRUE, n=2001, col="tomato3") > ## > (kk <- 100 * 2/1.24)# 161.29 [1] 161.2903 > curve(qgamma (x, 100, log=TRUE), -1.1*kk, -.95*kk, log='y', n=2001,col=2) > curve(qgammaAppr(x, 100, log=TRUE), add = TRUE, n=2001, col="tomato3") > abline(v = -kk, col='blue', lty=2)# exactly: kink is at a * 2 / 1.24 = a / .62 > curve(qgammaAppr(x - 100/.62, 100,log=TRUE), -1e-3, +1e-3) > > showProc.time() Time (user system elapsed): 0.16 0 0.16 > > p.qgammaLog <- function(alpha, xl.f = 1.5, xr.f = 0.4, n = 2001) + { + ## Purpose: + ## ---------------------------------------------------------------------- + ## Arguments: + ## ---------------------------------------------------------------------- + ## Author: Martin Maechler, Date: 30 Mar 2004, 18:44 + kk <- -alpha / .62 # = (alpha * 2) / (-1.24) + curve(qgamma(x, alpha, log=TRUE), xl.f*kk, xr.f*kk, log='y', + n=n, col=2, lwd=3.6, lty = 4, + main= paste("qgamma(x, alpha=",formatC(alpha,digits=10),", log = TRUE)")) + lines(kk, qgamma(kk, alpha, log=TRUE), type = 'h', lty = 3) + curve(qgamma (exp(x), alpha), add = TRUE, col="orange", n=n, lwd= 2) + curve(qgammaAppr(x, alpha, log=TRUE), add = TRUE, col=3, n=n,lwd = .4) + } > > showProc.time() Time (user system elapsed): 0.01 0 0.01 > > p.qgammaLog(25) > p.qgammaLog(16)# ~ [-25, -20] > p.qgammaLog(12, 1.2, 0.8)# small problem remaining > p.qgammaLog(11, 1.2, 0.8)# even smaller > p.qgammaLog(10.5, 1.1, 0.9)# even smaller > p.qgammaLog(10.25, 1.1, 0.9)# even smaller > ## 2019-08: __nothing__ visible from here on: > p.qgammaLog(10.18, 1.02, 0.98)# even smaller > p.qgammaLog(10.15, 1.02, 0.98)# even smaller > p.qgammaLog(10.14, 1.001, 0.999)# even smaller > p.qgammaLog(10.139, 1.0002, 0.9998)# > p.qgammaLog(10.138, 1.0002, 0.9998)# > p.qgammaLog(10.137, 1.00001, 0.99999)# > p.qgammaLog(10.13699, 1.0000001, 0.9999999)# > p.qgammaLog(10.1369899, 1.0000001, 0.9999999)# > p.qgammaLog(10.1369894, 1.0000001, 0.9999999)# > p.qgammaLog(10.1369893, 1.0000001, 0.9999999)# even smaller at -16.34998 > > showProc.time() Time (user system elapsed): 0.5 0 0.5 > > ##-- here is the boundary --- for 64-bit AMD Opteron --- > ## and for 32-bit AMD Athlon > > p.qgammaLog(10.1369892, 1.0000001, 0.9999999)# no more > p.qgammaLog(10.136989, 1.0000001, 0.9999999)# > p.qgammaLog(10.136988, 1.0000001, 0.9999999)# > p.qgammaLog(10.136985, 1.0000001, 0.9999999)# > p.qgammaLog(10.13698, 1.0000001, 0.9999999)# > p.qgammaLog(10.13697, 1.0000001, 0.9999999)# > p.qgammaLog(10.13695, 1.0000001, 0.9999999)# > p.qgammaLog(10.1368, 1.000001, 0.999999)# > p.qgammaLog(10.1365, 1.000001, 0.999999)# > p.qgammaLog(10.136, 1.000001, 0.999999)# > p.qgammaLog(10.125, 1.1, 0.9)# --- see it now > p.qgammaLog(10, 1.2, 0.8) > p.qgammaLog(9) > > showProc.time() Time (user system elapsed): 0.42 0 0.42 > > ## For large alpha: show difference to see problem better > ## ---> for alpha >= 10, the x problem starts *roughly* at x = -0.8*alpha > ## > > sfsmisc::mult.fig(2) > curve(qgammaAppr(x, 5, log=TRUE), - 8.1, -8, n=2001) > curve(qgammaAppr(x- 5/.62, 5, log=TRUE), -1e-15, 0) > > ## is the kink from pgamma() ? : no: this looks fine, > curve(pgamma(x, 1e5, log=TRUE), 1, 2e5, log='x', n=2001,col=2) > ## and this does too: > curve( dgamma(x, 1e5), .5e5, 2e5); par(new=TRUE) > curve( dgamma(x, 1e5, log=TRUE), .5e5, 2e5, col=2, yaxt="n") > axis(4,col.axis=2); par(new=TRUE) > curve( pgamma(x, 1e5), .5e5, 2e5, n=2001, col=3); par(new=TRUE) > curve( pgamma(x, 1e5, log=TRUE), .5e5, 2e5, n=2001, col=4); par(new=TRUE) > curve(-pgamma(x, 1e5, log=TRUE,lower=FALSE), .5e5, 2e5, n=2001, col=4) > ## all looking nice > > > x <- 10^seq(2,6, length=4001) > qx <- qgamma(pgamma(x, 1e5, log=TRUE), 1e5, log=TRUE) > plot(x, qx, type ='l', col=2, asp = 1); abline(0,1, lty=3) > > showProc.time() Time (user system elapsed): 0.07 0 0.07 > > ###------------- Approximations of qgamma() ------ > ## > > ## source("/u/maechler/R/MM/NUMERICS/dpq-functions/qchisqAppr.R") > ##--> qchisqAppr() > ##--> qchisqWH [ = Wilson Hilferty ] > ##--> qchisqKG [ = Kennedy & Gentle's improvements "a la AS 91" ] > ## dyn.load("/u/maechler/R/MM/NUMERICS/dpq-functions/qchisq_appr.so") > > ## Consider the two different implementations of > ## lgamma1p(a) := lgamma(1+a) == log(gamma(1+a) == log(a*gamma(a)) "stable": > > lseq <- if(requireNamespace("sfsmisc")) sfsmisc::lseq else + function(from, to, length) exp(seq(log(from), log(to), length.out = length)) > > if(require("Rmpfr")) { ##---------------- MPFR numbers ------------------------- + + .mpfr.all.eq <- Rmpfr::all.equal + AllEq <- function(target, current, ...) + .mpfr.all.eq(target, current, ..., + formatFUN = function(x, ...) Rmpfr::format(x, digits = 9)) + + print(gammaE <- Const("gamma",200)); pi. <- Const("pi",200) + print(a0 <- (gammaE^2 + pi.^2/6)/2) + print(psi2.1 <- -2*zeta(mpfr(3,200)))# == psigamma(1,2) =~ -2.4041138 + print(a1 <- (psi2.1 - gammaE*(pi.^2/2 + gammaE^2))/6) + + x <- lseq(1e-30, 0.8, length = if(doExtras) 1000 else 125) + x. <- mpfr(x, 200) + xct. <- log(x. * gamma(x.)) ## using MPFR arithmetic .. no overflow ... + xc2. <- log(x.) + lgamma(x.)## (ditto) + print(AllEq(xct., xc2., tol = 0)) # 3.15779......e-57 + xct <- as.numeric(xct.) + stopifnot(exprs = { + AllEq(xct., xc2., tol = 1e-45) + AllEq(xct , xc2., tol = 1e-15) + ## + all.equal(lgamma1p(x), lgamma1p(x, tol= 1e-16), tol=0) + ## -> no difference; i.e., default tol = 1e-14 seems fine enough! + }) + showProc.time() + + m.appr <- cbind(log(x*gamma(x)), lgamma(1+x), log(x) + lgamma(x), + lgamma1p.(x, k=1, cut=3e-6), + lgamma1p.(x, k=2, cut=1e-4), + lgamma1p.(x, k=3, cut=8e-4), + lgamma1p(x))#, tol= 1e-14), # = default + + eMat <- m.appr - xct # absolute error + ## Relative errors: + str(reMat. <- m.appr/xct. - 1) + str(reMat <- as(reMat., "array")) # as(., "matrix") fails in older versions + + matplot(x, eMat , log="x", type="l", lty=1) #-> problematic log(x) + lgamma(x) for "large" + matplot(x, abs( eMat), log="xy", type="l", lty=1) #-> but good for small; lgamma1p is much better + matplot(x, abs(reMat), log="xy", type="l", lty=1) + abline(v= 3.47548562941137e-08, col = "gray80", lwd=3)#<- the cutoff value of lgamma1p() + ##---> should use earlier cutoff! + ## zoom in: + + matplot(x, abs(reMat), log="xy", type="l", col=1:7, lty=1, + lwd=2, xlim=c(8e-9, 1e-3), ylim = c(1e-18, 1e-7), axes=FALSE, frame=TRUE, + main = expression(lgamma1p(x) == log(Gamma(x+1)) ~~~ "approximations" + ~~~ abs(rel.Err(.)))) + eaxis(1); eaxis(2) + abline(v= 3.47548562941137e-08, col = "gray80", lwd=3)#<- the cutoff value of lgamma1p() + abline(h= c(1,2,4)*.Machine$double.eps, lty=3, col="skyblue") + legend("topright", col=1:7, lty=1,lwd=2, + c("log(x*gamma(x))", "lgamma(1+x)", "log(x) + lgamma(x)", + "lgamma1p.(x, k=1, c=3e-6)", + "lgamma1p.(x, k=2, c=1e-4)", + "lgamma1p.(x, k=3, c=8e-4)", + "lgamma1p(x)"), bty="n", ncol=2) + abline(v = c(3e-6, 1e-4, 8e-4), col=4:6, lty=2, lwd=1/2) + + ## FIXME: do the same for the lgaamma1p_series() + + ## rm(x., xct., xc2., reMat., eMat, AllEq) + detach("package:Rmpfr") + showProc.time() + + } ## if( MPFR ) ---------------------------------------------------------------- Loading required package: Rmpfr Loading required package: gmp Attaching package: 'gmp' The following objects are masked from 'package:base': %*%, apply, crossprod, matrix, tcrossprod C code of R package 'Rmpfr': GMP using 64 bits per limb Attaching package: 'Rmpfr' The following object is masked from 'package:gmp': outer The following object is masked from 'package:DPQ': log1mexp The following objects are masked from 'package:stats': dbinom, dgamma, dnbinom, dnorm, dpois, dt, pnorm The following objects are masked from 'package:base': cbind, pmax, pmin, rbind 1 'mpfr' number of precision 200 bits [1] 0.57721566490153286060651209008240243104215933593992359880576723 1 'mpfr' number of precision 200 bits [1] 0.98905599532797255539539565150063470793918352072821409044319567 1 'mpfr' number of precision 200 bits [1] -2.404113806319188570799476323022899981529972584680997763584544 1 'mpfr' number of precision 200 bits [1] -0.90747907608088628901656016735627511492861144907256376094133062 [1] "Mean relative difference: 3.18105335e-57" Time (user system elapsed): 0.79 0.08 0.87 Class 'mpfrArray' [package "Rmpfr"] of dimension c(125, 7) and precision 200 -1 -1 ... num [1:125, 1:7] -1.00 -1.00 -1.00 3.64e+13 -1.00 ... - attr(*, "dimnames")=List of 2 ..$ : NULL ..$ : NULL Time (user system elapsed): 0.07 0 0.06 Warning messages: 1: In xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 348 y values <= 0 omitted from logarithmic plot 2: In xy.coords(x, y, xlabel, ylabel, log) : 1 y value <= 0 omitted from logarithmic plot > > > ## ../R/qchisqAppr.R -- talks about the "small shape" qgamma() approximation > ## ----------------- --> .qgammaApprBnd() : > curve(.qgammaApprBnd, 1e-18, 1e-15, col=2) > abline(h=0, col="gray70", lty=2) > eps.c <- .Machine$double.eps > axis(3,at=(1:3)* eps.c, + label=expression(epsilon[c], 2*epsilon[c], 3*epsilon[c])) > (rt.b <- uniroot(.qgammaApprBnd, c(1,3)*eps.c, tol=1e-12)) $root [1] 2.220446e-16 $f.root [1] 1.281676e-16 $iter [1] 0 $init.it [1] NA $estim.prec [1] 4.440892e-16 > rt.b$root ## 3.954775e-16 [1] 2.220446e-16 > rt.b$root / eps.c ## 1.781072 [1] 1 > ##==> for a < 1.781*eps, bnd > 0 ==> we have log(p) < bnd for all p > ## otherwise, we should effectively 'test' > curve(.qgammaApprBnd, 1e-16, 1e-10, log="x", col=2) > showProc.time() Time (user system elapsed): 0 0 0 > > > ## source ("/u/maechler/R/MM/NUMERICS/dpq-functions/beta-gamma-etc/qgamma-fn.R") > ## ##--> qchisqAppr.R() -- which has 'kind = ' argument! > ## ##--> qgamma.R() > > p.qchi.appr <- + function(x, qm= { m <- cbind(qchisq(x, df, log=TRUE), + sapply(knds, function(kind) + qchisqAppr.R(x,df,log=TRUE,kind=kind))) + colnames(m) <- c("True", "default", knds[-1]) + m }, + df, + knds = list(NULL,"chi.small", "WH", "p1WH", "df.small"), + call = match.call(), main = deparse(call), log = "y", ...) + { + ## Purpose: + ## ---------------------------------------------------------------------- + ## Arguments: + ## ---------------------------------------------------------------------- + ## Author: Martin Maechler, Date: 25 Mar 2004, 22:08 + + col <- c(2,1,3:6) + lty <- c(1,3,1,1,1,1) + lwd <- c(2,2,1,1,1,1) + matplot(x, qm, col=col, lty=lty, lwd=lwd, log = log, + type = 'l', main = main, ...) + y0 <- c( .02, .98) %*% par('usr')[3:4] + if(par("ylog")) y0 <- 10^y0 + legend(min(x), y0, colnames(qm), col=col, lty=lty, lwd=lwd) + invisible(list(x=x, qm=qm, call=match.call())) + } > > > pD.chi.appr <- function(pqr, err.kind=c("relative", "absolute"), + type = "l", log = "y", + lwds = c(2, rep(1, k-1)), + cols = seq(along=lwds), + ltys = rep(1,k), + ...) + { + ## Purpose: Plot Difference from "True" qchisq() + ## ---------------------------------------------------------------------- + ## Arguments: pqr: a list as resulting from p.chi.appr() + ## ---------------------------------------------------------------------- + ## Author: Martin Maechler, Date: 31 Mar 2004, 09:38 + err.kind <- match.arg(err.kind) + if(!is.list(pqr) || !is.numeric(k <- ncol(pqr$qm)-1) || k <= 0) + stop("invalid first argument 'pqr'") + with(pqr, { + err <- abs(if(err.kind == "relative") + (1- qm[,-1] / qm[,1]) else (qm[,-1] - qm[,1])) + matplot(x, err, ylab = "", + main = paste(err.kind,"Error from", deparse(call)), + type= type, log= log, lty=ltys, col=cols, lwd=lwds, ...) + legend(par("xaxp")[1], par("yaxp")[2], colnames(qm)[-1], + lty=ltys, col=cols, lwd=lwds) + }) + } > > ## if(FALSE) # you can manually set > ## do.pqchi <- TRUE # before source()ing this file > ## if(!exists("do.pqchi") || !is.logical(do.pqchi)) > ## do.pqchi <- interactive() > > ## if(do.pqchi) { #------------- FIXME look at speed up or cache indeed !?? <<<<< > > ## pqFile <- "/u/maechler/R/MM/NUMERICS/dpq-functions/pqchi.rda" > ## ## ls -l ................... 1325446 Nov 2 2009 pqchi.rda > ## if(file.exists(pqFile)) { > ## attach(pqFile) ## it loads more than we create here __FIXME__ > ## print(ls(2, all.names=TRUE)) > ## ## [1] "pq.1" "pq.25" "pq.25." "pq.31" "pq.33" "pq.33." "pq.33.2" > ## ## [8] "pq.33.3" "pq.33.4" "pq.5" "pq.5." "pq1" "pq1." "pq1.95" > ## ## [15] "pq1.95." "pq1.95.2" "pq10" "pq10." "pq10.2" "pq100" "pq2" > ## ## [22] "pq2." "pq2.05" "pq2.05." "pq2.05.2" "pq2.5" "pq2.5." "pq2.5.2" > ## ## [29] "pq200" "pq2L" "pq4" "pq4." "pq4.2" "pqFile" > ## } > > showProc.time() Time (user system elapsed): 0 0 0 > x <- seq(-300, -.01, length=501) > > all.equal(qchisqAppr.R(x, 200, log=TRUE), + qchisqAppr (x, 200, log=TRUE),tol=0) [1] "Mean relative difference: 1.968318e-16" > ## 4.48 e-16 / TRUE (Opteron) > > all.equal(qchisqAppr.R(x, 2, log=TRUE), + qchisqAppr (x, 2, log=TRUE),tol=0) [1] "Mean relative difference: 1.535551e-16" > ## 3.90 e-16 / TRUE (Opteron) > > all.equal(qchisqAppr.R(x, 0.1, log=TRUE), + qchisqAppr (x, 0.1, log=TRUE),tol=0) [1] TRUE > ## 7.15 e-15 / 1.179e-8 !!!!! (Opteron) > > pq200 <- p.qchi.appr(x = seq(-300, -.01, length=501), df = 200) Warning message: In xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 939 y values <= 0 omitted from logarithmic plot > pq100 <- p.qchi.appr(x = seq(-160, -.01, length=501), df = 100) Warning message: In xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 941 y values <= 0 omitted from logarithmic plot > ## after (slow) computing, quickly repeat: > with(pq200, p.qchi.appr(x=x, qm=qm, call=call)) Warning message: In xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 939 y values <= 0 omitted from logarithmic plot > with(pq100, p.qchi.appr(x=x, qm=qm, call=call)) Warning message: In xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 941 y values <= 0 omitted from logarithmic plot > > ## this "hangs forever" -- before I introduced 'maxit' (for 'nu.small'): > pq10 <- p.qchi.appr(x = seq(-12, -.005, length=501), df = 10) Warning message: In xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 891 y values <= 0 omitted from logarithmic plot > ## want to see the jump: > pq10. <- p.qchi.appr(x = seq(-10, -4, length=501), df = 10) Warning message: In xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 1002 y values <= 0 omitted from logarithmic plot > pq10.2<- p.qchi.appr(x = seq(-8.5,-7.5, length=501), df = 10) Warning message: In xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 1002 y values <= 0 omitted from logarithmic plot > with(pq10.2, p.qchi.appr(x=x, qm=qm, call=call)) Warning message: In xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 1002 y values <= 0 omitted from logarithmic plot > > > pq2.5 <- p.qchi.appr(x = seq(-3.4, -.01, length=501), df = 2.5) Warning message: In xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 244 y values <= 0 omitted from logarithmic plot > pq2.5. <- p.qchi.appr(x = seq(-2.1, -1.8, length=901), df = 2.5)#the jump Warning message: In xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 901 y values <= 0 omitted from logarithmic plot > ## what about p1WH (which is fantastic for df=2)? > pq2.5.2<- p.qchi.appr(x = seq(-0.5, -1e-3, length=901), df = 2.5) > with(pq2.5, p.qchi.appr(x=x, qm=qm, call=call)) Warning message: In xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 244 y values <= 0 omitted from logarithmic plot > with(pq2.5.2, p.qchi.appr(x=x, qm=qm, call=call)) > pD.chi.appr(pq2.5.2)# nothing special > > pq2.05 <- p.qchi.appr(x = seq(-3.4, -.01, length=501), df = 2.05) Warning message: In xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 81 y values <= 0 omitted from logarithmic plot > pq2.05. <- p.qchi.appr(x = seq(-2.5, -1.5, length=901), df = 2.05)#the jump > ## ^^ the jump from chi.small to WH is much too late here > ## what about p1WH (which is fantastic for df=2)? > pq2.05.2<- p.qchi.appr(x = seq(-0.4, -1e-5, length=901), df = 2.05) > pD.chi.appr(pq2.05.2) # p1WH is starting to become better > # and the jump (WH -> p1WH) is too late > > with(pq2.05, p.qchi.appr(x=x, qm=qm, call=call)) Warning message: In xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 81 y values <= 0 omitted from logarithmic plot > with(pq2.05.2, p.qchi.appr(x=x, qm=qm, call=call)) > > > ## Here, all are 'ok' (but "nu.small"): > pq2L <- p.qchi.appr(seq(-300, -.01, length=201), df = 2) There were 50 or more warnings (use warnings() to see the first 50) > > pq2 <- p.qchi.appr(x = seq(-5, -.001, length=501), df = 2) > pq2. <- p.qchi.appr(x = seq(-2.5, -1, length=901), df = 2) > with(pq2., p.qchi.appr(x=x, qm=qm, call=call)) > > pq4 <- p.qchi.appr(x = seq(-8, -0.01, length = 501), df = 4) Warning message: In xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 845 y values <= 0 omitted from logarithmic plot > summary(warnings()) 1 identical warnings: In xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 845 y values <= 0 omitted from logarithmic plot > pq4.2 <- p.qchi.appr(x = seq(-0.1, -1e-04, length = 901), df = 4) > > pq1.95 <- p.qchi.appr(x = seq(-3., -.01, length=501), df = 1.95) > pq1.95. <- p.qchi.appr(x = seq(-2.2, -1.5, length=901), df = 1.95)#the jump -1.57 > ## ^^ the jump from chi.small to WH is *much* too late here > ## what about p1WH (which is fantastic for df=2)? > pq1.95.2<- p.qchi.appr(x = seq(-0.4, -1e-7, length=901), df = 1.95) > pD.chi.appr(pq1.95.2) # p1WH is starting to become better > # and the jump (WH -> p1WH) is too late > > with(pq1.95, p.qchi.appr(x=x, qm=qm, call=call)) > with(pq1.95.2, p.qchi.appr(x=x, qm=qm, call=call)) > > > pq1 <- p.qchi.appr(x = seq(-4, -.001, length=501), df = 1) There were 50 or more warnings (use warnings() to see the first 50) > pq1. <- p.qchi.appr(x = seq(-1.8, -.5, length=901), df = 1) > with(pq1., p.qchi.appr(x=x, qm=qm, call=call)) > > pq.5 <- p.qchi.appr(x = seq(-1.5, -.001, length=501), df = 0.5) Warning message: In xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 243 y values <= 0 omitted from logarithmic plot > pq.5. <- p.qchi.appr(x = seq(-0.8, -0.2, length=901), df = 0.5, ylim=c(.04,.6)) Warning message: In xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 40 y values <= 0 omitted from logarithmic plot > > pq.33 <- p.qchi.appr(x= seq(-0.9, -.001,length=501), df= 0.33) Warning message: In xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 246 y values <= 0 omitted from logarithmic plot > pq.33. <- p.qchi.appr(x= seq(-0.4, -0.02,length=901), df= 0.33) > pq.33.2<- p.qchi.appr(x= seq(-0.4, -0.2, length=901), df= 0.33, ylim=c(.15,.60)) > with(pq.33.2, p.qchi.appr(x=x, qm=qm, call=call,ylim=c(.15,.45))) > > pq.33.3<- p.qchi.appr(x= seq(-0.4, -0.005, length=901), df= 0.33, ylim=c(.15, 4.00)) > with(pq.33.3, p.qchi.appr(x=x, qm=qm, call=call))#,ylim=c(.15, 8))) > > pq.33.4<- p.qchi.appr(x= seq(-0.003, -1e-6, length=901), df= 0.33,ylim=c(5,25)) > with(pq.33.4, p.qchi.appr(x=x, qm=qm, call=call,ylim=c(5,25))) > > ## nu <= 0.32 is the "magic" border of Best & Roberts > > pq.31 <- p.qchi.appr(x = seq(-0.45,-.010, length=501), df = 0.31) Warning message: In xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 27 y values <= 0 omitted from logarithmic plot > with(pq.31, p.qchi.appr(x=x, qm=qm, call=call)) Warning message: In xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 27 y values <= 0 omitted from logarithmic plot > > pq.25 <- p.qchi.appr(x = seq(-0.3, -0.02, length=901), df = 0.25) > > pq.1 <- p.qchi.appr(x = seq(-0.16, -0.01, length=901), df = 0.1) Warning message: In xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 226 y values <= 0 omitted from logarithmic plot > with(pq.1, p.qchi.appr(x=x, qm=qm, call=call)) Warning message: In xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 226 y values <= 0 omitted from logarithmic plot > showProc.time() Time (user system elapsed): 3.48 0.06 3.55 > > ## if(!file.exists(pqFile)) # don't overwrite for now (as it contains pq2L , > ## save(list=ls(pat="^pq"), file = pqFile) > > ##}## end if(do.pqchi){ only if interactive } ====================================== > > pD.chi.appr(pq2L, "abs") Warning messages: 1: In xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 363 y values <= 0 omitted from logarithmic plot 2: In xy.coords(x, y, xlabel, ylabel, log) : 180 y values <= 0 omitted from logarithmic plot > pD.chi.appr(pq2L, "rel") Warning messages: 1: In xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 363 y values <= 0 omitted from logarithmic plot 2: In xy.coords(x, y, xlabel, ylabel, log) : 180 y values <= 0 omitted from logarithmic plot > ## --> want only much smaller x-range: > pD.chi.appr(pq2,"abs")#--> fantastic p1WH Warning messages: 1: In xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 340 y values <= 0 omitted from logarithmic plot 2: In xy.coords(x, y, xlabel, ylabel, log) : 1 y value <= 0 omitted from logarithmic plot > pD.chi.appr(pq2) # (ditto) Warning messages: 1: In xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 340 y values <= 0 omitted from logarithmic plot 2: In xy.coords(x, y, xlabel, ylabel, log) : 1 y value <= 0 omitted from logarithmic plot > > pD.chi.appr(pq4.2)# p1WH: only at very right > pD.chi.appr(pq4.2, xlim=c(-.016,0))# p1WH: only at very right > > > ## no Newton step here, eg: > (qgA100 <- qgammaAppr(1e-100, 100)) [1] 3.799269 > (qg.100 <- qgamma (1e-100, 100)) [1] 3.950799 > all.equal(qgA100, qg.100) [1] "Mean relative difference: 0.03988396" > ## too much different > dgamma(1e-100, 100, log = TRUE)# -23154.7 i.e., "non-log" is 0 [1] -23154.73 > > qgamma.R(1e-100, 100, verbose = TRUE)#-> final Newton fails! qgamma(p= 1e-100, alpha= 100, scale= 1, l.t.= 1, log.p= 0): (p1,df)=(-230.259,200) ==> small chi-sq., ch0 = 7.59854 ==> ch = 7.59854: Ph.II iter; ch=7.59854, p2=9.76738e-101 it=2, ch=2.45916e+09, p2=-1 --> Del became NA it=1: p=-230.259, x = 3.79927, p.=-234.019; p1:=D{p}=-3.76094 new t= 3.94774005, p.= -230.332933 it=2, d{p}=-0.074424 new t= 3.95079758, p.= -230.258539 it=3, d{p}=-2.99766e-05 new t= 3.95079881, p.= -230.258509 it=4, d{p}=-4.88853e-12 new t= 3.95079881, p.= -230.258509 it=5, d{p}=0 [1] 3.950799 > > ## but here, the final Newton iterations do work : > x <- qgamma.R(log(1e-100), 100, log = TRUE, verbose = TRUE) qgamma(p=-230.259, alpha= 100, scale= 1, l.t.= 1, log.p= 1): (p1,df)=(-230.259,200) ==> small chi-sq., ch0 = 7.59854 it=1: p=-230.259, x = 3.79927, p.=-234.019; p1:=D{p}=-3.76094 new t= 3.94774005, p.= -230.332933 it=2, d{p}=-0.074424 new t= 3.95079758, p.= -230.258539 it=3, d{p}=-2.99766e-05 new t= 3.95079881, p.= -230.258509 it=4, d{p}=-4.88853e-12 new t= 3.95079881, p.= -230.258509 it=5, d{p}=0 > pgamma(x, 100) # = 1e-100 ! perfect [1] 1e-100 > showProc.time() Time (user system elapsed): 0.06 0 0.06 > > ###--> Use this to devise an improved final Newton algorithm !!! > > > ## From: Prof Brian Ripley > ## To: skylab.gupta@gmail.com > ## Cc: R-bugs@biostat.ku.dk, r-devel@stat.math.ethz.ch > ## Subject: Re: [Rd] qgamma inaccuracy (PR#12324) > ## Date: Tue, 12 Aug 2008 20:50:50 +0100 (BST) > > ## This is a really extreme usage. AFAICS the code works well enough down to > ## shape=1e-10 or so, e.g. > > qgamma(1e-10, 5e-11, lower.tail=FALSE) [1] 0.08237203 > ## [1] 0.08237203 > ## in R 2.9.. .. 2.10.0 -- with an accuracy warning {which is *wrong*!} > > ## I would be interested to know what substantive problem you were trying to > ## solve here that required such values. > > ## I am pretty sure that a completely different algorithm will be required. > > ## MM: It looks like this is basis for a new algo: > a <- 1e-14 > gammaE <- 0.57721566490153286060651209008240243079 # Euler's gamma > curve(pgamma(x, a, lower.tail=FALSE)/a + log(x) + gammaE, 1e-300, 1e-1, log="",n=1000) > curve(pgamma(x, a, lower.tail=FALSE)/a + log(x) + gammaE - x, 1e-300, 1e-1, log="",n=1000) > ## ==> Q = 1 - P = pgamma(x,a, lower=FALSE) ~= a*(-log(x) - gammaE + x - x^2/4) > ## i.e., Q ~= -a(log(x) + gammaE { -x + x^2/4 } > ## -Q/a - gammaE ~= log(x) { -x + x^2/4 } > ## ==> x ~= exp(- (Q/a + gammaE)) > ## e.g., example below: > Q <- 1e-100; a <- 5e-101 > ## MM: Find inverse : > str(r.a <- uniroot(function(x) pgamma(x,a,lower.tail=FALSE) - Q, + int = c(0.01, 0.1), tol=1e-20)) List of 5 $ root : num 0.0824 $ f.root : num 0 $ iter : int 8 $ init.it : int NA $ estim.prec: num 5.95e-11 > dput(x0 <- r.a$root) ## 0.0823720296207203 0.0823720296207203 > (x1 <- exp(- (Q/a + gammaE)))## 0.07598528 .. not so good [1] 0.07598528 > qgammaApprSmallP(Q, a, lower.tail=FALSE)## ~= 0.07598528 -- the same!! [1] 0.07598528 > > pgamma(x0, a, lower.tail=FALSE) ## 1.00000e-100. [1] 1e-100 > pgamma(x1, a, lower.tail=FALSE) ## 1.03728e-100 ... hmm "close" [1] 1.037283e-100 > ## > > ## MM: -- now look at the bigger picture > p.qg.2a <- function(l2x.min= -15, l2x.max = -100, n = 501, + do.offset = FALSE, + type = "o", log = "x", cex = 0.6, ...) { + x.log <- any("x" == strsplit(log,"")[[1]]) + x <- if(x.log) 2^seq(l2x.min, l2x.max, length=n) + else seq(2^l2x.min,2^l2x.max, length=n) + if(do.offset) + plot(x, qgamma(2*x, x, lower.tail=FALSE) - 0.0823720296206873, + type=type, cex=cex, log=log, ...) + else plot(x, qgamma(2*x, x, lower.tail=FALSE), + type=type, cex=cex, log=log, ...) + } > p.qg.2a() # was "very bad" in R <= 2.10.0, now --> 0.082372... "perfect" smooth > ## still a little ---but acceptable--- remaining inaccuracy ...zooming in: > p.qg.2a(-43,-55, do.offset=TRUE) > p.qg.2a(,-1024) > p.qg.2a(,-1024, log="", pch=".")## linear in x !! > ## zoom in at the limit > p.qg.2a(-30,-1024, do.offset=TRUE, ylim = 1e-11*c(-1,1)) > p.qg.2a(-33,-1024, do.offset=TRUE, ylim = 1e-12*c(-1,1)) > p.qg.2a(-33,-1024, do.offset=TRUE, ylim = 1e-13*c(-1,1)) > > a <- 2^-(7:900) > qg <- qgamma(2*a, a, lower.tail=FALSE) > re <- 1-pgamma(qg, a, lower.tail=FALSE)/(2*a) > plot(a, re, log="x", type="b", col=2) > stopifnot(abs(re) < 2e-12) # but really, *should be a bit better > > showProc.time() Time (user system elapsed): 0.14 0.03 0.17 > ## For completeness we may write that in due course, but for now (R 2.7.2) I > ## suggest just issuing a warning for miniscule 'shape'. > > ## On Thu, 7 Aug 2008, skylab.gupta@gmail.com wrote: > > ## > Full_Name: > ## > Version: 2.7.1 (2008-06-23) > ## > OS: windows vista > ## > Submission from: (NULL) (216.82.144.137) > ## > > ## > > ## > Hello, > ## > > ## > I have been working with various probability distributions in R, and it seems > ## > the gamma distribution is inaccurate for some inputs. > > ## > For example, qgamma(1e-100, 5e-101, lower.tail=FALSE) gives: 1.0. However, it > ## > seems this is incorrect; I think the correct answer should be > ## > 0.082372029620717283. When I check these numbers using pgamma, I get: > > (qg <- qgamma(1e-100, 5e-101, lower.tail=FALSE))# 1 (wrong, originally) [1] 0.08237203 > ## 0.08237203 now (2009-11-04), i.e. ok > pgamma(qg, 5e-101, lower.tail=FALSE)# now -> 1e-100 : ok [1] 1e-100 > > pgamma(0.082372029620717283,5e-101, lower.tail=FALSE) [1] 1e-100 > ## 1.0000000000000166e-100. > > RE.pqgamma <- function(p, shape, lower.tail = TRUE, log.p = FALSE) { + ## Relative Error of pgamma(qgamma(*), ..): + 1 - pgamma(qgamma(p, shape, lower.tail=lower.tail, log.p=log.p), + shape=shape, lower.tail=lower.tail, log.p=log.p) / p + } > RE.qpgamma <- function(q, shape, lower.tail = TRUE, log.p = FALSE) { + ## Relative Error of qgamma(pgamma(*), ..): + 1 - qgamma(pgamma(q, shape, lower.tail=lower.tail, log.p=log.p), + shape=shape, lower.tail=lower.tail, log.p=log.p) / q + } > > ## Ok, how extreme can we get -- let a := alpha := shape --> 0 : > x <- 1e-100 > a <- 2^-(7:300)# is still "ok": > plot(a, (re <- RE.pqgamma(x, a, lower.tail=FALSE)), type="b", col=2, log="x")# oops! > a <- 2^-(7:400) > plot(a, (re <- RE.pqgamma(x, a, lower.tail=FALSE)), type="b", col=2, log="x")# oops! > ## Oops! > ## but, it is clear > qgamma(x, 2^-400, lower.tail=FALSE)## is exactly 0 [1] 0 > > ## -> it goes to 0 quickly . .. zooming in: > curve(qgamma(1e-100, x, lower.tail=FALSE), 1e-110, 1e-70, log="xy", col=2, axes=FALSE) Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 18 y values <= 0 omitted from logarithmic plot > eaxis(1);eaxis(2) > > ## from when on is it exactly 0: > uniroot(function(u) qgamma(1e-100, 2^u,lower.tail=FALSE)-1e-315, c(-400, -300))$root [1] -341.6941 > ## -341.6941 > ## => use > a <- 2^-(7:341) > plot(a, (re <- RE.pqgamma(x, a, lower.tail=FALSE)), + type="l", col=2, log="x")# small glips > ## zoom in: > curve(abs(RE.pqgamma(1e-100, x, lower.tail=FALSE)), 2^-341, 1e-90, log="xy", n=2000) Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 7 y values <= 0 omitted from logarithmic plot > curve(RE.pqgamma(1e-100, x, lower.tail=FALSE), 1e-100, 10e-100, n=2000) > curve(RE.pqgamma(1e-100, x, lower.tail=FALSE), 4e-100, 6e-100, n=2000) > ## Ok: at least here is a problem > RE.pqgamma(1e-100, 5e-100, lower.tail=FALSE)# -0.1538 [1] 3.874678e-14 > > ## more general > curve(RE.pqgamma(x, 5*x, lower.tail=FALSE), 1e-100, 1e-20, n=10000, log="x") > ## problem *everywhere* , starting quite early: (lesser problem at ~ 1e-16 !) > curve(RE.pqgamma(x, 5*x, lower.tail=FALSE), 1e-25, 1e-15, n=1000, log="x") > curve(RE.pqgamma(x, 5*x, lower.tail=FALSE), 1e-21, 10e-21,n=1000) > curve(RE.pqgamma(x, 5*x, lower.tail=FALSE), 2e-21, 6e-21, n=1000) > curve(RE.pqgamma(x, 5*x, lower.tail=FALSE), 4e-21, 4.5e-21, n=1000) > ## and indeed, it's qgamma() that jumps here > curve(qgamma(x, 5*x, lower.tail=FALSE), 4e-21, 4.5e-21, ylim=c(.97, 1.1)) > > ## well, looking at pgamma(), finally reveals the buglet is there first: > ## There's a jump at x = 1 !! > curve(pgamma(x, 1e-30, lower=FALSE), .9999, 1.0001) > curve(pgamma(x, 1e-20, lower=FALSE), .9999, 1.0001) > curve(pgamma(x, 1e-17, lower=FALSE), .9999, 1.0001) > curve(pgamma(x, 1e-15, lower=FALSE), .9999, 1.0001) > curve(pgamma(x, 1e-13, lower=FALSE), .9999, 1.0001) > curve(pgamma(x, 1e-12, lower=FALSE), .9999, 1.0001)# still clearly visible > curve(pgamma(x, 1e-11, lower=FALSE), .9999, 1.0001)# barely visible > ## for larger alpha == shape, must zoom in more and more: > curve(pgamma(x, 1e-11, lower=FALSE), .99999, 1.00001)# > curve(pgamma(x, 1e-10, lower=FALSE), .999999, 1.000001) > curve(pgamma(x, 1e-9, lower=FALSE), .9999999, 1.0000001) > curve(pgamma(x, 1e-8, lower=FALSE), .99999999, 1.00000001) > curve(pgamma(x, 1e-7, lower=FALSE), .999999999, 1.000000001) > curve(pgamma(x, 1e-6, lower=FALSE), .9999999999, 1.0000000001) > curve(pgamma(x, 1e-5, lower=FALSE), .99999999999, 1.00000000001) > curve(pgamma(x, 1e-4, lower=FALSE), .999999999999, 1.000000000001) > ## now we get close to noise level: > curve(pgamma(x, 1e-3, lower=FALSE), .9999999999999, 1.0000000000001) > curve(pgamma(x, 1e-2, lower=FALSE), .99999999999999, 1.00000000000001) > > showProc.time() Time (user system elapsed): 0.24 0.02 0.25 > > del.pgamma <- function(a, eps = 1e-13) + { + ## Purpose: *relative* jump size at x = 1 of pgamma(x, a, lower=FALSE) + ## ---------------------------------------------------------------------- + ## Author: Martin Maechler, Date: 5 Nov 2009, 16:08 + stopifnot(a > 0, length(a) == 1, eps > 0, length(eps) == 1) + pp <- pgamma(1+c(-eps,eps), a, lower.tail = FALSE) + (pp[2] - pp[1])*2/(pp[2] + pp[1]) + } > > a <- lseq(1e-300, 1e-3, length=400) > dpa <- sapply(a, del.pgamma) > plot(a, -dpa, log="xy", type="l", col=2, yaxt="n");eaxis(2) > ## ok, it remains constant all the way to 1e-300 > ## --> focus > a <- lseq(1e-40, 1e-5, length=400) > dpa <- sapply(a, del.pgamma) > plot(a, -dpa, log="xy", type="l", col=2, axes=FALSE) > eaxis(1, at = 10^-seq(5,40, by=5)) > eaxis(2) > > > xm <- .Machine$double.xmin > pgamma(xm, shape= 1e-20)# is "practically 1" --> *most* qgamma() will be exactly 0 [1] 1 > ## how close to 1 ? ---> use upper_tail [possibly on log scale:] > pgamma(xm, shape= 1e-20, lower.tail=FALSE) # 7.078192e-18 [1] 7.078192e-18 > pgamma(xm, shape= 1e-20, lower.tail=FALSE, log=TRUE) # -39.48951 [1] -39.48951 > > ## Where is the 'boundary' (from which on qgamma() must return 0, since it can't give > ## xm = 2.2e-308 ) ? > a <- 2^-(7:1030) > plot(a, (p <- pgamma(xm, a, lower.tail=FALSE, log=TRUE)), + cex=.5,type ="l", col=2, log="x") > summary(lm. <- lm(p ~ log(a) + a + I(a^2))) ## coeff. of log(a) *is* 1 Call: lm(formula = p ~ log(a) + a + I(a^2)) Residuals: Min 1Q Median 3Q Max -0.0081566 -0.0000098 0.0000171 0.0000440 0.0107065 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 6.562e+00 3.327e-05 197240 <2e-16 *** log(a) 1.000e+00 8.024e-08 12463221 <2e-16 *** a -3.403e+02 2.039e-01 -1669 <2e-16 *** I(a^2) 1.551e+04 2.911e+01 533 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.0005225 on 1020 degrees of freedom Multiple R-squared: 1, Adjusted R-squared: 1 F-statistic: 5.249e+13 on 3 and 1020 DF, p-value: < 2.2e-16 > summary(lm. <- lm(p ~ offset(log(a)) + a + I(a^2))) Call: lm(formula = p ~ offset(log(a)) + a + I(a^2)) Residuals: Min 1Q Median 3Q Max -0.0081788 0.0000179 0.0000179 0.0000179 0.0107351 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 6.562e+00 1.639e-05 400476.3 <2e-16 *** a -3.404e+02 2.033e-01 -1674.4 <2e-16 *** I(a^2) 1.552e+04 2.907e+01 533.7 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.0005232 on 1021 degrees of freedom Multiple R-squared: 0.9999, Adjusted R-squared: 0.9999 F-statistic: 8.426e+06 on 2 and 1021 DF, p-value: < 2.2e-16 > > p.l <- pgamma(xm, a, lower.tail=FALSE, log=TRUE) - log(a) > dput(mean(tail(p.l))) ## 6.56218869790132 6.56218869790132 > ##=> pgamma(xm, a) ~= log(a) + 6.5621886979 + > ## ok, to get this better, now need different a: > al <- seq(1e-300, 1e-3, length=200) > plot(al, (pl <- pgamma(xm, al, lower.tail=FALSE, log=TRUE)) - (log(al)+6.56218869790132), + cex=.5,type ="l", col=2) > summary(lm.. <- lm(pl ~ offset(log(al) + 6.56218869790132) + 0 + al + I(al^2))) Call: lm(formula = pl ~ offset(log(al) + 6.56218869790132) + 0 + al + I(al^2)) Residuals: Min 1Q Median 3Q Max -1.201e-05 -4.268e-06 -1.336e-06 3.075e-06 5.247e-06 Coefficients: Estimate Std. Error t value Pr(>|t|) al -3.539e+02 2.032e-03 -174123 <2e-16 *** I(al^2) 2.075e+04 2.617e+00 7929 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 4.153e-06 on 198 degrees of freedom Multiple R-squared: 1, Adjusted R-squared: 1 F-statistic: 2.217e+11 on 2 and 198 DF, p-value: < 2.2e-16 > confint(lm..) 2.5 % 97.5 % al -353.8628 -353.8547 I(al^2) 20745.6597 20755.9814 > coef(lm..) al I(al^2) -353.8587 20750.8205 > ## al I(al^2) > ## -353.8587 20750.8205 > ##=> pgamma(xm, a) ~= log(a) + 6.5621886979 - 353.858745 * a > > ## > Similarly, for example: -- now (2009-11-04) ok > qgamma(1e-100,0.005,lower.tail=FALSE) # = 109.36757177 instead of 219.5937.. [1] 219.5937 > pgamma(109.36757177007101, 0.005, lower.tail=FALSE)# = 1.4787306506694758e-52. [1] 1.478731e-52 > > ## > This looks completely wrong. The correct value, I think, should be > ## > 219.59373661415756. In fact, > pgamma(219.59373661415756, 0.005, lower.tail=FALSE)# = 9.9999999999999558e-101. [1] 1e-100 > ## > > ## > In fact, when I do the following in R, the results are completely wrong, > ## > > a <- 5*10^-(1:40) > z1 <- qgamma(1e-100,a,lower.tail=FALSE) > (y <- pgamma(z1,a,lower.tail=FALSE)) [1] 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 [11] 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 [21] 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 [31] 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 > ## The value of y that I get should be close to 1e-100, but they are not: > ## [1] 1.000000e-100 1.871683e-51 1.478731e-52 1.444034e-53 1.440606e-54 > ## [6] 1.440264e-55 1.440230e-56 1.440226e-57 1.440226e-58 1.440226e-59 > summary(abs(1 - y/1e-100)) Min. 1st Qu. Median Mean 3rd Qu. Max. 2.220e-16 4.663e-15 1.221e-14 1.488e-14 1.940e-14 5.262e-14 > plot(a, abs(1 - y/1e-100), log="xy", type="b") > stopifnot(abs(1 - y/1e-100) < 2e-13)# max (32b Linux, P.M) = 4.186e-14 > > ## > The correct values of z1 should be: > z1true <- c(226.97154111939946, 222.15218724493326, 219.59373661415756, + 217.27485383840451, 214.98015408183574, 212.68797118872064, + 210.39614286838227, 208.10445550564617, 205.81289009100664, + 203.52144711679352) > all.equal(z1[1:10], z1true, tol=0)# 1.307e-16 on 32-bit (Pentium M); 1.686e-16 (64b Lnx, 2019) [1] "Mean relative difference: 1.761977e-16" > stopifnot(all.equal(z1[1:10], z1true, tol = 1e-15)) > showProc.time() Time (user system elapsed): 0.06 0 0.07 > ##> > ##> With these values of z1true, we get the expected values: > > (y <- pgamma(z1true, a, lower.tail=FALSE)) [1] 1.000000e-100 1.000000e-100 1.000000e-100 1.000000e-100 1.000000e-100 [6] 1.000000e-100 1.000000e-100 1.000000e-100 1.000000e-100 1.000000e-100 [11] 5.869617e-112 7.428625e-111 9.705940e-111 9.970232e-111 9.997024e-111 [16] 9.999703e-111 9.999970e-111 9.999997e-111 1.000000e-110 1.000000e-110 [21] 5.869617e-122 7.428625e-121 9.705940e-121 9.970232e-121 9.997024e-121 [26] 9.999703e-121 9.999970e-121 9.999997e-121 1.000000e-120 1.000000e-120 [31] 5.869617e-132 7.428625e-131 9.705940e-131 9.970232e-131 9.997024e-131 [36] 9.999703e-131 9.999970e-131 9.999997e-131 1.000000e-130 1.000000e-130 > ## [1] 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 > > ## > I am using the precompiled binary version of R, under Windows Vista. > ## > ----------- > ## >> version > ## > _ > ## > platform i386-pc-mingw32 > ## > arch i386 > ## > os mingw32 > ## > system i386, mingw32 > ## > status > ## > major 2 > ## > minor 7.1 > ## > year 2008 > ## > month 06 > ## > day 23 > ## > svn rev 45970 > ## > language R > ## > version.string R version 2.7.1 (2008-06-23) > ## > ------------ > ## > > ## > So, it seems qgamma is inaccurate for small probability values in the upper > ## > tail, when the shape parameter is also small. > > > ###_-- MM: Still wrong: > > (xm <- 2^-1074.9999) # is less than .Machine $ double.xmin == the really x > 0 [1] 4.940656e-324 > > pgamma(xm, .00001)# 0.992589 [1] 0.992589 > qgamma(.99, .00001)##--> NaN -- should give 0 or "xmin" or so [1] 0 > ## FIXME -- ok, now > > ## but > curve(qgamma(x, .001, lower=FALSE), .4, .8, n=1001, log="y") Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 688 y values <= 0 omitted from logarithmic plot > curve(qgamma(x, 1e-5, lower=FALSE), .002, .2, n=1001, log="xy") Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 716 y values <= 0 omitted from logarithmic plot > curve(qgamma(x, 1e-7, lower=FALSE), 1e-5, .04, n=1001, log="xy") Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 759 y values <= 0 omitted from logarithmic plot > curve(qgamma(x, 1e-12, lower=FALSE), 1e-12, 1e-2, n=1001, log="xy") Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 713 y values <= 0 omitted from logarithmic plot > > ## or > curve(qgamma(x, 1e-121, lower=FALSE), 7e-119, 8e-119, + n=2001, log="y", yaxt="n") Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 1117 y values <= 0 omitted from logarithmic plot > try( # reveals eaxis() bug ? -- for the *subnormal* numbers + eaxis(2, at = 10^-seq(304,324, by=2)) + ) Error in eaxis(2, at = 10^-seq(304, 324, by = 2)) : invalid 'log=TRUE' for at <= 0: not a true log scale plot? > > curve(qgamma(x, .001, lower=FALSE), .4, .6, n=1001, log="y") Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 376 y values <= 0 omitted from logarithmic plot > curve(qgamma(x, .001, lower=FALSE), .5, .55, n=1001, log="y") Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 503 y values <= 0 omitted from logarithmic plot > ## gives a *warning* from axis() because of subnormal y-range {was error, fixed in R-devel (2018-08)} > curve(qgamma(x, .001, lower=FALSE), .52, .53, n=1001, log="y") Warning messages: 1: In xy.coords(x, y, xlabel, ylabel, log) : 514 y values <= 0 omitted from logarithmic plot 2: In plot.window(...) : R_pretty(): very small range 'cell'=2.5e-308, corrected to 2.50321e-308 3: In plot.window(...) : axis(2, *): range of values ( 0) is small wrt |M| = 1e-307 --> not pretty() 4: In axis(side = side, at = at, labels = labels, ...) : CreateAtVector [log-axis()]: small axp[0] = 4.94066e-324 > > curve(qgamma(x, .001, lower=FALSE)*1e100, .522, .526, n=1001, log="y") Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 285 y values <= 0 omitted from logarithmic plot > > showProc.time() Time (user system elapsed): 0.03 0 0.03 > > proc.time() user system elapsed 11.65 0.53 12.17