* using log directory 'd:/Rcompile/CRANpkg/local/3.5/DPQ.Rcheck' * using R version 3.5.3 (2019-03-11) * using platform: x86_64-w64-mingw32 (64-bit) * using session charset: ISO8859-1 * checking for file 'DPQ/DESCRIPTION' ... OK * this is package 'DPQ' version '0.3-5' * package encoding: UTF-8 * checking package namespace information ... OK * checking package dependencies ... OK * checking if this is a source package ... OK * checking if there is a namespace ... OK * checking for hidden files and directories ... OK * checking for portable file names ... OK * checking whether package 'DPQ' can be installed ... WARNING Found the following significant warnings: Warning: unable to re-encode 'pnbeta.R' lines 12, 35 Warning: unable to re-encode 'pnchisq.R' lines 263, 264, 265, 266, 648, 652, 654, 667, 675, 676, 716 Warning: unable to re-encode 'qnchisq.R' lines 135, 136, 146, 147, 148, 149, 243, 248, 250 See 'd:/Rcompile/CRANpkg/local/3.5/DPQ.Rcheck/00install.out' for details. * checking installed package size ... OK * checking package directory ... OK * checking 'build' directory ... OK * checking DESCRIPTION meta-information ... OK * checking top-level files ... OK * checking for left-over files ... OK * checking index information ... OK * checking package subdirectories ... OK * checking R files for non-ASCII characters ... OK * checking R files for syntax errors ... OK * loading checks for arch 'i386' ** checking whether the package can be loaded ... OK ** checking whether the package can be loaded with stated dependencies ... OK ** checking whether the package can be unloaded cleanly ... OK ** checking whether the namespace can be loaded with stated dependencies ... OK ** checking whether the namespace can be unloaded cleanly ... OK ** checking loading without being on the library search path ... OK ** checking use of S3 registration ... OK * loading checks for arch 'x64' ** checking whether the package can be loaded ... OK ** checking whether the package can be loaded with stated dependencies ... OK ** checking whether the package can be unloaded cleanly ... OK ** checking whether the namespace can be loaded with stated dependencies ... OK ** checking whether the namespace can be unloaded cleanly ... OK ** checking loading without being on the library search path ... OK ** checking use of S3 registration ... OK * checking dependencies in R code ... OK * checking S3 generic/method consistency ... OK * checking replacement functions ... OK * checking foreign function calls ... OK * checking R code for possible problems ... [14s] OK * checking Rd files ... OK * checking Rd metadata ... OK * checking Rd cross-references ... OK * checking for missing documentation entries ... OK * checking for code/documentation mismatches ... OK * checking Rd \usage sections ... OK * checking Rd contents ... OK * checking for unstated dependencies in examples ... OK * checking line endings in C/C++/Fortran sources/headers ... OK * checking pragmas in C/C++ headers and code ... OK * checking compiled code ... OK * checking sizes of PDF files under 'inst/doc' ... OK * checking installed files from 'inst/doc' ... OK * checking files in 'vignettes' ... OK * checking examples ... ** running examples for arch 'i386' ... [8s] OK ** running examples for arch 'x64' ... [7s] OK * checking for unstated dependencies in 'tests' ... OK * checking tests ... ** running tests for arch 'i386' ... [124s] ERROR Running 'chisq-nonc-ex.R' [31s] Running 'dnchisq-tst.R' [1s] Running 'pnbeta-tst.R' [1s] Running 'pnt-precision-problem.R' [22s] Running 'ppois-ex.R' [3s] Running 'qPoisBinom-ex.R' [1s] Running 'qbeta-dist.R' [18s] Running 'qbeta-tst.R' [1s] Running 'qgamma-ex.R' [24s] Running 't-nonc-tst.R' [9s] Running 'wienergerm-accuracy.R' [9s] Running 'wienergerm-pchisq-tst.R' [1s] Running 'wienergerm_nchisq.R' [3s] Running the tests in 'tests/qgamma-ex.R' failed. Complete output: > 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="Matrix", "test-tools-1.R", mustWork=TRUE)) Loading required package: tools > ##--> showProc.time(), assertError(), relErrV(), ... > showProc.time() Time elapsed: 1.83 0.08 1.97 > > (doExtras <- DPQ:::doExtras()) [1] FALSE > (sdir <- system.file("safe", package="DPQ")) ## save directory (to read from) [1] "D:/temp/RtmpcJZjAr/RLIBS_dec02f1a3c6f/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 elapsed: 0.79 0.01 0.83 > > require(graphics) > if(!dev.interactive(orNone=TRUE)) pdf("qgamma-appr.pdf") > eaxis <- sfsmisc::eaxis > > showProc.time() Time elapsed: 0.09 0 0.09 > ## 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 6.46592989151e-08 7.2 [2,] 0.000605300028912 0.02041606754775 -1.99820160418e-11 10.7 [3,] 0.001012316063255 0.01855615147677 -2.59145555105e-04 3.6 [4,] 0.001248285290785 0.01838201076117 3.94217658517e-10 9.4 [5,] 0.001682388899865 0.01720736646288 5.53088974600e-04 3.3 [6,] 0.001746787400790 0.01731189518997 -5.86897856980e-08 7.2 [7,] 0.002664451237518 0.01599317398629 1.48421342013e-04 3.8 [8,] 0.002664451237518 0.01618024201222 7.75564700239e-08 7.1 [9,] 0.003159421860255 0.01557612780310 -7.92117005632e-06 5.1 [10,] 0.003159421860255 0.01568183691729 -4.52237560733e-08 7.3 [11,] 0.004055462418244 0.01493858731306 -7.17404703909e-06 5.1 [12,] 0.004400694140827 0.01459101672970 9.07907026435e-04 3.0 [13,] 0.004458811277768 0.01457506850867 9.03139988525e-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.63408184153e-10 9.6 [19,] 0.011642124851282 0.01201471267173 -3.00271327153e-04 3.5 [20,] 0.014753716559535 0.01155624353203 1.52961088240e-10 9.8 [21,] 0.015499213434879 0.01125420134457 1.00492447767e-04 4.0 [22,] 0.015499213434879 0.01135920381800 -9.55739012376e-08 7.0 [23,] 0.018603016576955 0.01071716109330 -2.07537936111e-03 2.7 [24,] 0.018603016576955 0.01073655493589 2.14388784341e-04 3.7 [25,] 0.022624242394389 0.01033379525113 -3.37865668776e-09 8.5 [26,] 0.022624242394389 0.01034206121729 4.43079739565e-08 7.4 [27,] 0.023730217356634 0.01016252135853 -1.07732682664e-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 -1.78633419168e-09 8.7 [31,] 0.081818424963746 0.00689856335721 -2.30915286892e-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 -6.63889809793e-09 8.2 [36,] 0.110869751789691 0.00603268830251 8.95578944338e-10 9.0 [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 8.96455243371e-10 9.0 [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.79807193976e-10 9.3 [52,] 0.190094017358772 0.00450373552308 1.78565421871e-06 5.7 [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.58750790291e-11 10.8 [56,] 0.201518808589718 0.00439976887947 -9.97180116258e-11 10.0 [57,] 0.210803673024037 0.00427351441034 -1.70232716812e-10 9.8 [58,] 0.213058614771766 0.00426179831847 -1.39808165045e-11 10.9 [59,] 0.214780951412088 0.00419869272965 -1.91879370171e-08 7.7 [60,] 0.232805106603566 0.00395399315002 -9.17581020055e-08 7.0 [61,] 0.249102914025652 0.00380019404026 -1.15818687974e-10 9.9 [62,] 0.249102914025652 0.00382493512126 -1.39672717836e-11 10.9 [63,] 0.252076511947811 0.00374903834738 2.14569900181e-07 6.7 [64,] 0.253082914021191 0.00375259362798 3.65436092498e-09 8.4 [65,] 0.253922058700076 0.00371237348323 -5.09014937222e-06 5.3 [66,] 0.254289278570932 0.00374343873151 -1.05664899053e-09 9.0 [67,] 0.260017499519858 0.00366179605930 -2.90430199001e-07 6.5 [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 -6.07109695849e-10 9.2 [72,] 0.284601233201101 0.00344936339590 -2.63026489478e-10 9.6 [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 -4.48927431229e-07 6.3 [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.82253101499e-11 10.7 [81,] 0.303182283998467 0.00328723648952 -2.23643326080e-11 10.7 [82,] 0.322319846303892 0.00306134512927 -1.15130830540e-05 4.9 [83,] 0.322319846303892 0.00310689001755 8.57751647487e-11 10.1 [84,] 0.325071272052651 0.00302343293053 -2.47088704493e-04 3.6 [85,] 0.325071272052651 0.00304146419577 -7.71376615361e-06 5.1 [86,] 0.331888412404218 0.00300837121343 -4.96098895297e-09 8.3 [87,] 0.362278153188527 0.00278204202032 4.53939330569e-10 9.3 [88,] 0.385389476781711 0.00260981704384 -1.77430203863e-09 8.8 [89,] 0.425333956955001 0.00232995789362 1.82823025607e-08 7.7 [90,] 0.439503709203564 0.00222452690840 -4.53585193982e-06 5.3 [91,] 0.439503709203564 0.00224964327069 -3.02331937263e-10 9.5 [92,] 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 6.465930e-08 7.2 [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.868979e-08 7.2 [5,] 0.0026644512 0.015993174 1.484213e-04 3.8 [6,] 0.0026644512 0.016180242 7.755647e-08 7.1 [7,] 0.0031594219 0.015576128 -7.921170e-06 5.1 [8,] 0.0031594219 0.015681837 -4.522376e-08 7.3 [9,] 0.0040554624 0.014938587 -7.174047e-06 5.1 [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 -3.002713e-04 3.5 [17,] 0.0154992134 0.011254201 1.004924e-04 4.0 [18,] 0.0154992134 0.011359204 -9.557390e-08 7.0 [19,] 0.0186030166 0.010717161 -2.075379e-03 2.7 [20,] 0.0186030166 0.010736555 2.143888e-04 3.7 [21,] 0.0226242424 0.010333795 -3.378657e-09 8.5 [22,] 0.0226242424 0.010342061 4.430797e-08 7.4 [23,] 0.0237302174 0.010162521 -1.077327e-06 6.0 [24,] 0.0447535254 0.008396264 1.222242e-05 4.9 [25,] 0.0818184250 0.006860077 -1.786334e-09 8.7 [26,] 0.0828003091 0.006812347 4.179976e-09 8.4 [27,] 0.0908216581 0.006552698 -7.160336e-09 8.1 [28,] 0.1022947605 0.006235631 -6.638898e-09 8.2 [29,] 0.1274058577 0.005623691 6.605415e-09 8.2 [30,] 0.1352296342 0.005400734 -2.347626e-05 4.6 [31,] 0.1377322800 0.005330921 2.992858e-04 3.5 [32,] 0.1381129175 0.005351387 -2.053358e-06 5.7 [33,] 0.1411006360 0.005273058 4.315938e-05 4.4 [34,] 0.1429052994 0.005236800 3.481808e-04 3.5 [35,] 0.1456245572 0.005269240 -1.945018e-09 8.7 [36,] 0.1546068729 0.005068069 -4.599247e-07 6.3 [37,] 0.1546068729 0.005073662 2.723010e-07 6.6 [38,] 0.1697410365 0.004841818 5.314010e-09 8.3 [39,] 0.1773275763 0.004659561 5.534044e-05 4.3 [40,] 0.1900940174 0.004503736 1.785654e-06 5.7 [41,] 0.1901476415 0.004534687 5.662356e-09 8.2 [42,] 0.2147809514 0.004198693 -1.918794e-08 7.7 [43,] 0.2328051066 0.003953993 -9.175810e-08 7.0 [44,] 0.2520765119 0.003749038 2.145699e-07 6.7 [45,] 0.2530829140 0.003752594 3.654361e-09 8.4 [46,] 0.2539220587 0.003712373 -5.090149e-06 5.3 [47,] 0.2542892786 0.003743439 -1.056649e-09 9.0 [48,] 0.2600174995 0.003661796 -2.904302e-07 6.5 [49,] 0.2703239068 0.003519992 -1.561648e-04 3.8 [50,] 0.2716993561 0.003550681 5.130930e-09 8.3 [51,] 0.2755161961 0.003468040 -4.351715e-04 3.4 [52,] 0.2901885431 0.003366135 -5.644431e-08 7.2 [53,] 0.2905790220 0.003344235 1.026676e-07 7.0 [54,] 0.2905790220 0.003367649 2.260616e-08 7.6 [55,] 0.2918501987 0.003335528 -1.273388e-06 5.9 [56,] 0.2965211365 0.003303089 -4.489274e-07 6.3 [57,] 0.2980341749 0.003304623 8.424704e-09 8.1 [58,] 0.3005567833 0.003239225 4.660033e-05 4.3 [59,] 0.3223198463 0.003061345 -1.151308e-05 4.9 [60,] 0.3250712721 0.003023433 -2.470887e-04 3.6 [61,] 0.3250712721 0.003041464 -7.713766e-06 5.1 [62,] 0.3318884124 0.003008371 -4.960989e-09 8.3 [63,] 0.3853894768 0.002609817 -1.774302e-09 8.8 [64,] 0.4253339570 0.002329958 1.828230e-08 7.7 [65,] 0.4395037092 0.002224527 -4.535852e-06 5.3 [66,] 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.912e-04 -1.441e-04 -1.969e-05 7.710e-05 1.082e-03 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 6.184e-06 1.133e-05 0.546 0.587 log(p) -2.725e-03 3.660e-06 -744.530 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.0002557 on 90 degrees of freedom Multiple R-squared: 0.9998, Adjusted R-squared: 0.9998 F-statistic: 5.543e+05 on 1 and 90 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.094448581017 -7.463565388e-08 7.1 2 9.685654347e-16 0.094212475025 1.225666370e-06 5.9 3 1.056228096e-15 0.093976369034 4.405535633e-07 6.4 4 1.151824906e-15 0.093740263042 -3.124345540e-07 6.5 5 1.256073967e-15 0.093504157050 1.062787617e-06 6.0 6 1.369758374e-15 0.093268051059 3.686723420e-07 6.4 7 1.493732098e-15 0.093031945067 -2.933092969e-07 6.5 8 1.628926404e-15 0.092795839075 1.156863211e-06 5.9 9 1.776356839e-15 0.092559733084 -1.520860305e-06 5.8 10 1.937130869e-15 0.092323627092 -1.717405684e-08 7.8 11 2.112456192e-15 0.092087521100 1.507978934e-06 5.8 12 2.303649813e-15 0.091851415109 -1.062629703e-06 6.0 13 2.512147934e-15 0.091615309117 5.160568806e-07 6.3 14 2.739516748e-15 0.091379203125 6.832590693e-08 7.2 15 2.987464197e-15 0.091143097134 -3.472479155e-07 6.5 16 3.257852808e-15 0.090906991142 1.306469104e-06 5.9 17 3.552713679e-15 0.090670885150 -1.081912813e-06 6.0 18 3.874261739e-15 0.090434779159 6.253705814e-07 6.2 19 4.224912384e-15 0.090198673167 3.330732011e-07 6.5 20 4.607299625e-15 0.089962567175 7.294750648e-08 7.1 21 5.024295868e-15 0.089726461183 -1.550036943e-07 6.8 22 5.479033495e-15 0.089490355192 -3.507775794e-07 6.5 23 5.974928394e-15 0.089254249200 1.485185219e-06 5.8 24 6.515705616e-15 0.089018143208 -6.457821511e-07 6.2 25 7.105427358e-15 0.088782037217 1.243792663e-06 5.9 26 7.748523477e-15 0.088545931225 -8.120437036e-07 6.1 27 8.449824769e-15 0.088309825233 1.131156029e-06 5.9 28 9.214599250e-15 0.088073719242 -8.495397843e-07 6.1 29 1.004859174e-14 0.087837613250 1.147297759e-06 5.9 30 1.095806699e-14 0.087601507258 -7.582479513e-07 6.1 31 1.194985679e-14 0.087365401267 1.292240281e-06 5.9 32 1.303141123e-14 0.087129295275 -5.381458088e-07 6.3 33 1.421085472e-14 0.086893189283 -3.797838852e-07 6.4 34 1.549704695e-14 0.086657083292 -1.892109929e-07 6.7 35 1.689964954e-14 0.086420977300 3.357565381e-08 7.5 36 1.842919850e-14 0.086184871308 2.885788447e-07 6.5 37 2.009718347e-14 0.085948765317 -1.348411111e-06 5.9 38 2.191613398e-14 0.085712659325 8.952460000e-07 6.0 39 2.389971358e-14 0.085476553333 -6.665524503e-07 6.2 40 2.606282246e-14 0.085240447341 -2.772836227e-07 6.6 41 2.842170943e-14 0.085004341350 1.442152372e-07 6.8 42 3.099409391e-14 0.084768235358 -1.299326189e-06 5.9 43 3.379929908e-14 0.084532129366 1.083914191e-06 6.0 44 3.685839700e-14 0.084296023375 -2.844134916e-07 6.5 45 4.019436694e-14 0.084059917383 2.714025346e-07 6.6 46 4.383226796e-14 0.083823811391 8.594620383e-07 6.1 47 4.779942715e-14 0.083587705400 -3.905788137e-07 6.4 48 5.212564492e-14 0.083351599408 2.673437002e-07 6.6 49 5.684341886e-14 0.083115493416 9.575175537e-07 6.0 50 6.198818782e-14 0.082879387425 -1.742195652e-07 6.8 51 6.759859815e-14 0.082643281433 -1.262888024e-06 5.9 52 7.371679400e-14 0.082407175441 1.378141781e-06 5.9 53 8.038873388e-14 0.082171069450 3.647252201e-07 6.4 54 8.766453592e-14 0.081934963458 -6.056185433e-07 6.2 55 9.559885430e-14 0.081698857466 2.942142103e-07 6.5 56 1.042512898e-13 0.081462751475 1.226316300e-06 5.9 57 1.136868377e-13 0.081226645483 3.743139556e-07 6.4 58 1.239763756e-13 0.080990539491 -4.346107259e-07 6.4 59 1.351971963e-13 0.080754433499 -1.200456962e-06 5.9 60 1.474335880e-13 0.080518327508 -1.231656597e-07 6.9 61 1.607774678e-13 0.080282221516 9.864073824e-07 6.0 62 1.753290718e-13 0.080046115524 3.389273247e-07 6.5 63 1.911977086e-13 0.079810009533 -2.654694915e-07 6.6 64 2.085025797e-13 0.079573903541 -8.267823173e-07 6.1 65 2.273736754e-13 0.079337797549 4.280209399e-07 6.4 66 2.479527513e-13 0.079101691558 -5.255530455e-08 7.3 67 2.703943926e-13 0.078865585566 1.272196181e-06 5.9 68 2.948671760e-13 0.078629479574 8.723627990e-07 6.1 69 3.215549355e-13 0.078393373583 5.156188723e-07 6.3 70 3.506581437e-13 0.078157267591 2.019651378e-07 6.7 71 3.823954172e-13 0.077921161599 -6.859768598e-08 7.2 72 4.170051594e-13 0.077685055608 -2.960688799e-07 6.5 73 4.547473509e-13 0.077448949616 -4.804477280e-07 6.3 74 4.959055026e-13 0.077212843624 -6.217335187e-07 6.2 75 5.407887852e-13 0.076976737633 -7.199255436e-07 6.1 76 5.897343520e-13 0.076740631641 -7.750230990e-07 6.1 77 6.431098711e-13 0.076504525649 -7.870254888e-07 6.1 78 7.013162874e-13 0.076268419657 -7.559320068e-07 6.1 79 7.647908344e-13 0.076032313666 -6.817419649e-07 6.2 80 8.340103188e-13 0.075796207674 -5.644546737e-07 6.2 81 9.094947018e-13 0.075560101682 -4.040694535e-07 6.4 82 9.918110051e-13 0.075323995691 -2.005856121e-07 6.7 83 1.081577570e-12 0.075087889699 4.599752257e-08 7.3 84 1.179468704e-12 0.074851783707 3.356806249e-07 6.5 85 1.286219742e-12 0.074615677716 6.684643596e-07 6.2 86 1.402632575e-12 0.074379571724 1.044349404e-06 6.0 87 1.529581669e-12 0.074143465732 -1.904324292e-07 6.7 88 1.668020638e-12 0.073907359741 2.770713794e-07 6.6 89 1.818989404e-12 0.073671253749 -8.551932531e-07 6.1 90 1.983622010e-12 0.073435147757 -2.960718908e-07 6.5 91 2.163155141e-12 0.073199041766 3.061522486e-07 6.5 92 2.358937408e-12 0.072962935774 9.514798076e-07 6.0 93 2.572439484e-12 0.072726829782 1.875341071e-08 7.7 94 2.805265149e-12 0.072490723791 -8.599822618e-07 6.1 95 3.059163338e-12 0.072254617799 -7.452604778e-08 7.1 96 3.336041275e-12 0.072018511807 7.540343130e-07 6.1 97 3.637978807e-12 0.071782405815 2.630667650e-08 7.6 98 3.967244020e-12 0.071546299824 -6.474386178e-07 6.2 99 4.326310282e-12 0.071310193832 3.212443331e-07 6.5 100 4.717874816e-12 0.071074087840 -2.500146477e-07 6.6 101 5.144878969e-12 0.070837981849 -7.672983209e-07 6.1 102 5.610530299e-12 0.070601875857 3.415009425e-07 6.5 103 6.118326675e-12 0.070365769865 -7.330663121e-08 7.1 104 6.672082550e-12 0.070129663874 -4.341459592e-07 6.4 105 7.275957614e-12 0.069893557882 8.147632093e-07 6.1 106 7.934488041e-12 0.069657451890 -9.939312102e-07 6.0 107 8.652620563e-12 0.069421345899 3.519772401e-07 6.5 108 9.435749632e-12 0.069185239907 2.015231014e-07 6.7 109 1.028975794e-11 0.068949133915 1.050243966e-07 7.0 110 1.122106060e-11 0.068713027924 6.247828421e-08 7.2 111 1.223665335e-11 0.068476921932 7.388190371e-08 7.1 112 1.334416510e-11 0.068240815940 1.392323935e-07 6.9 113 1.455191523e-11 0.068004709949 2.585268818e-07 6.6 114 1.586897608e-11 0.067768603957 -1.074910069e-06 6.0 115 1.730524113e-11 0.067532497965 6.589363856e-07 6.2 116 1.887149926e-11 0.067296391974 -5.557286085e-07 6.3 117 2.057951587e-11 0.067060285982 -2.152397722e-07 6.7 118 2.244212120e-11 0.066824179990 1.791772732e-07 6.7 119 2.447330670e-11 0.066588073998 -8.518538748e-07 6.1 120 2.668833020e-11 0.066351968007 -3.441463585e-07 6.5 121 2.910383046e-11 0.066115862015 2.174792778e-07 6.7 122 3.173795216e-11 0.065879756023 -6.299664652e-07 6.2 123 3.461048225e-11 0.065643650032 4.492504890e-08 7.3 124 3.774299853e-11 0.065407544040 7.737245498e-07 6.1 125 4.115903175e-11 0.065171438048 1.098156234e-07 7.0 126 4.488424239e-11 0.064935332057 -4.892612773e-07 6.3 127 4.894661340e-11 0.064699226065 -1.023513942e-06 6.0 128 5.337666040e-11 0.064463120073 -6.281881859e-08 7.2 129 5.820766091e-11 0.064227014082 9.517654395e-07 6.0 130 6.347590433e-11 0.063990908090 6.009614978e-07 6.2 131 6.922096451e-11 0.063754802098 3.149542980e-07 6.5 132 7.548599706e-11 0.063518696107 9.373602483e-08 7.0 133 8.231806350e-11 0.063282590115 -6.270114161e-08 7.2 134 8.976848478e-11 0.063046484123 -1.543650285e-07 6.8 135 9.789322680e-11 0.062810378132 -1.812634687e-07 6.7 136 1.067533208e-10 0.062574272140 -1.434043004e-07 6.8 137 1.164153218e-10 0.062338166148 -4.079537130e-08 7.4 138 1.269518087e-10 0.062102060156 1.265554705e-07 6.9 139 1.384419290e-10 0.061865954165 3.586403706e-07 6.4 140 1.509719941e-10 0.061629848173 6.554514619e-07 6.2 141 1.646361270e-10 0.061393742181 1.016980880e-06 6.0 142 1.795369696e-10 0.061157636190 8.972756516e-08 7.0 143 1.957864536e-10 0.060921530198 -7.618360451e-07 6.1 144 2.135066416e-10 0.060685424206 -1.952735760e-07 6.7 145 2.328306437e-10 0.060449318215 4.359719799e-07 6.4 146 2.539036173e-10 0.060213212223 -1.996387393e-07 6.7 147 2.768838580e-10 0.059977106231 -7.596117781e-07 6.1 148 3.019439882e-10 0.059741000240 7.654519618e-08 7.1 149 3.292722540e-10 0.059504894248 9.773499157e-07 6.0 150 3.590739391e-10 0.059268788256 -6.763557745e-07 6.2 151 3.915729072e-10 0.059032682265 3.646181882e-07 6.4 152 4.270132832e-10 0.058796576273 1.716164724e-07 6.8 153 4.656612873e-10 0.058560470281 5.417191207e-08 7.3 154 5.078072346e-10 0.058324364290 1.227028412e-08 7.9 155 5.537677160e-10 0.058088258298 4.589736990e-08 7.3 156 6.038879765e-10 0.057852152306 1.550389357e-07 6.8 157 6.585445080e-10 0.057616046314 3.396807524e-07 6.5 158 7.181478783e-10 0.057379940323 -6.657672476e-07 6.2 159 7.831458144e-10 0.057143834331 9.354081640e-07 6.0 160 8.540265665e-10 0.056907728339 9.185895744e-08 7.0 161 9.313225746e-10 0.056671622348 -6.652309521e-07 6.2 162 1.015614469e-09 0.056435516356 -9.234395604e-08 7.0 163 1.107535432e-09 0.056199410364 5.559531132e-07 6.3 164 1.207775953e-09 0.055963304373 4.705635348e-08 7.3 165 1.317089016e-09 0.055727198381 -3.754633768e-07 6.4 166 1.436295757e-09 0.055491092389 -7.116280665e-07 6.1 167 1.566291629e-09 0.055254986398 2.545663202e-07 6.6 168 1.708053133e-09 0.055018880406 8.553017516e-08 7.1 169 1.862645149e-09 0.054782774414 2.785924047e-09 8.6 170 2.031228938e-09 0.054546668423 6.311568379e-09 8.2 171 2.215070864e-09 0.054310562431 9.608511653e-08 7.0 172 2.415551906e-09 0.054074456439 2.720845634e-07 6.6 173 2.634178032e-09 0.053838350448 -6.486866118e-07 6.2 174 2.872591513e-09 0.053602244456 -2.948026030e-07 6.5 175 3.132583258e-09 0.053366138464 1.452384917e-07 6.8 176 3.416106266e-09 0.053130032472 6.714146586e-07 6.2 177 3.725290298e-09 0.052893926481 1.227572537e-07 6.9 178 4.062457877e-09 0.052657820489 -3.287850758e-07 6.5 179 4.430141728e-09 0.052421714497 4.666330337e-07 6.3 180 4.831103812e-09 0.052185608506 2.036972987e-07 6.7 181 5.268356064e-09 0.051949502514 3.778649948e-08 7.4 182 5.745183026e-09 0.051713396522 -3.113050129e-08 7.5 183 6.265166516e-09 0.051477290531 -3.084833278e-09 8.5 184 6.832212532e-09 0.051241184539 1.218923642e-07 6.9 185 7.450580597e-09 0.051005078547 3.437699647e-07 6.5 186 8.124915754e-09 0.050768972556 -4.487175920e-07 6.3 187 8.860283457e-09 0.050532866564 -2.762538487e-08 7.6 188 9.662207623e-09 0.050296760572 4.902705625e-07 6.3 189 1.053671213e-08 0.050060654581 1.026316276e-08 8.0 190 1.149036605e-08 0.049824548589 -3.619620201e-07 6.4 191 1.253033303e-08 0.049588442597 -6.264466410e-07 6.2 192 1.366442506e-08 0.049352336606 -7.832323377e-07 6.1 193 1.490116119e-08 0.049116230614 2.401634375e-07 6.6 194 1.624983151e-08 0.048880124622 -7.738734984e-07 6.1 195 1.772056691e-08 0.048644018630 -6.078122046e-07 6.2 196 1.932441525e-08 0.048407912639 -3.342184687e-07 6.5 197 2.107342426e-08 0.048171806647 4.686610289e-08 7.3 198 2.298073210e-08 0.047935700655 5.353999194e-07 6.3 199 2.506066606e-08 0.047699594664 9.197667672e-08 7.0 200 2.732885013e-08 0.047463488672 -2.330247564e-07 6.6 201 2.980232239e-08 0.047227382680 5.886113930e-07 6.2 202 3.249966302e-08 0.046991276689 -5.279761615e-07 6.3 203 3.544113383e-08 0.046755170697 5.191576196e-07 6.3 204 3.864883049e-08 0.046519064705 -3.498820569e-07 6.5 205 4.214684851e-08 0.046282958714 -8.357652814e-08 7.1 206 4.596146421e-08 0.046046852722 3.008301066e-07 6.5 207 5.012133212e-08 0.045810746730 -1.917574097e-07 6.7 208 5.465770025e-08 0.045574640739 -5.552768079e-07 6.3 209 5.960464478e-08 0.045338534747 1.941360427e-07 6.7 210 6.499932603e-08 0.045102428755 8.300436694e-08 7.1 211 7.088226765e-08 0.044866322764 1.007455417e-07 7.0 212 7.729766099e-08 0.044630216772 -7.200064360e-07 6.1 213 8.429369702e-08 0.044394110780 -4.391848569e-07 6.4 214 9.192292842e-08 0.044158004788 -2.969477242e-08 7.5 215 1.002426642e-07 0.043921898797 5.083973372e-07 6.3 216 1.093154005e-07 0.043685792805 2.298753619e-07 6.6 217 1.192092896e-07 0.043449686813 9.089335073e-08 7.0 218 1.299986521e-07 0.043213580822 9.137020085e-08 7.0 219 1.417645353e-07 0.042977474830 2.312248708e-07 6.6 220 1.545953220e-07 0.042741368838 -4.125440396e-07 6.4 221 1.685873940e-07 0.042505262847 1.135778693e-08 7.9 222 1.838458568e-07 0.042269156855 5.743894688e-07 6.2 223 2.004853285e-07 0.042033050863 3.701760187e-07 6.4 224 2.186308010e-07 0.041796944872 3.160064802e-07 6.5 225 2.384185791e-07 0.041560838880 4.117840361e-07 6.4 226 2.599973041e-07 0.041324732888 6.574119313e-07 6.2 227 2.835290706e-07 0.041088626897 1.687303717e-07 6.8 228 3.091906439e-07 0.040852520905 -1.591887169e-07 6.8 229 3.371747881e-07 0.040616414913 5.464674205e-07 6.3 230 3.676917137e-07 0.040380308922 -3.331954996e-07 6.5 231 4.009706570e-07 0.040144202930 -1.795109488e-07 6.7 232 4.372616020e-07 0.039908096938 1.344805564e-07 6.9 233 4.768371582e-07 0.039671990946 -2.420257981e-07 6.6 234 5.199946082e-07 0.039435884955 -4.473476847e-07 6.3 235 5.670581412e-07 0.039199778963 -4.816173824e-07 6.3 236 6.183812879e-07 0.038963672971 -3.449670454e-07 6.5 237 6.743495762e-07 0.038727566980 -3.752869437e-08 7.4 238 7.353834273e-07 0.038491460988 4.405657880e-07 6.4 239 8.019413140e-07 0.038255354996 -5.454495733e-07 6.3 240 8.745232040e-07 0.038019249005 2.846233194e-07 6.5 241 9.536743164e-07 0.037783143013 -3.274450626e-07 6.5 242 1.039989216e-06 0.037547037021 5.340015219e-08 7.3 243 1.134116282e-06 0.037310931030 -1.797782818e-07 6.7 244 1.236762576e-06 0.037074825038 5.685058231e-07 6.2 245 1.348699152e-06 0.036838719046 -7.036316796e-08 7.2 246 1.470766855e-06 0.036602613055 -5.064464588e-07 6.3 247 1.603882628e-06 0.036366507063 3.281364691e-08 7.5 248 1.749046408e-06 0.036130401071 -3.854574882e-09 8.4 249 1.907348633e-06 0.035894295080 1.616862507e-07 6.8 250 2.079978433e-06 0.035658189088 -2.268208801e-07 6.6 251 2.268232565e-06 0.035422083096 -4.023808344e-07 6.4 252 2.473525152e-06 0.035185977105 -3.652116916e-07 6.4 253 2.697398305e-06 0.034949871113 -1.155312188e-07 6.9 254 2.941533709e-06 0.034713765121 3.464431350e-07 6.5 255 3.207765256e-06 0.034477659129 -4.359180803e-07 6.4 256 3.498092816e-06 0.034241553138 4.610668849e-07 6.3 257 3.814697266e-06 0.034005447146 1.355166542e-07 6.9 258 4.159956866e-06 0.033769341154 4.357489369e-08 7.4 259 4.536465130e-06 0.033533235163 1.849733936e-07 6.7 260 4.947050303e-06 0.033297129171 -1.408841688e-07 6.9 261 5.394796609e-06 0.033061023179 -2.228060141e-07 6.7 262 5.883067419e-06 0.032824917188 -6.108725370e-08 7.2 263 6.415530512e-06 0.032588811196 3.439775413e-07 6.5 264 6.996185632e-06 0.032352705204 3.140542086e-07 6.5 265 7.629394531e-06 0.032116599213 -1.344452274e-07 6.9 266 8.319913732e-06 0.031880493221 -3.182555455e-07 6.5 267 9.072930260e-06 0.031644387229 4.235536309e-07 6.4 268 9.894100606e-06 0.031408281238 1.067806606e-07 7.0 269 1.078959322e-05 0.031172175246 6.478579662e-08 7.2 270 1.176613484e-05 0.030936069254 2.971861416e-07 6.5 271 1.283106102e-05 0.030699963263 -4.743406139e-07 6.3 272 1.399237126e-05 0.030463857271 3.167996901e-07 6.5 273 1.525878906e-05 0.030227751279 1.255002516e-07 6.9 274 1.663982746e-05 0.029991645287 2.292534416e-07 6.6 275 1.814586052e-05 0.029755539296 1.095587332e-08 8.0 276 1.978820121e-05 0.029519433304 9.794742695e-08 7.0 277 2.157918644e-05 0.029283327312 -1.157428633e-07 6.9 278 2.353226967e-05 0.029047221321 -1.397550675e-08 7.9 279 2.566212205e-05 0.028811115329 4.027337354e-07 6.4 280 2.798474253e-05 0.028575009337 -4.365435546e-08 7.4 281 3.051757812e-05 0.028338903346 -1.538802443e-07 6.8 282 3.327965493e-05 0.028102797354 7.146577696e-08 7.1 283 3.629172104e-05 0.027866691362 5.978360873e-08 7.2 284 3.957640242e-05 0.027630585371 3.936098777e-07 6.4 285 4.315837288e-05 0.027394479379 -4.939450760e-08 7.3 286 4.706453935e-05 0.027158373387 -1.258946654e-07 6.9 287 5.132424410e-05 0.026922267396 -3.862692297e-07 6.4 288 5.596948506e-05 0.026686161404 -2.704251889e-07 6.6 289 6.103515625e-05 0.026450055412 2.208878911e-07 6.7 290 6.655930986e-05 0.026213949421 2.101241836e-08 7.7 291 7.258344208e-05 0.025977843429 2.173189355e-07 6.7 292 7.915280485e-05 0.025741737437 -2.345979131e-07 6.6 293 8.631674575e-05 0.025505631445 -2.697597914e-07 6.6 294 9.412907870e-05 0.025269525454 1.109045821e-07 7.0 295 1.026484882e-04 0.025033419462 -1.036342423e-07 7.0 296 1.119389701e-04 0.024797313470 1.180052477e-07 6.9 297 1.220703125e-04 0.024561207479 -2.129809857e-07 6.7 298 1.331186197e-04 0.024325101487 -8.761589876e-08 7.1 299 1.451668842e-04 0.024088995495 1.024239626e-08 8.0 300 1.583056097e-04 0.023852889504 9.613828733e-08 7.0 301 1.726334915e-04 0.023616783512 1.855632931e-07 6.7 302 1.882581574e-04 0.023380677520 2.939531960e-07 6.5 303 2.052969764e-04 0.023144571529 -2.374739316e-08 7.6 304 2.238779402e-04 0.022908465537 1.742085209e-07 6.8 305 2.441406250e-04 0.022672359545 -1.221569534e-08 7.9 306 2.662372394e-04 0.022436253554 -1.074149700e-07 7.0 307 2.903337683e-04 0.022200147562 -9.628875608e-08 7.0 308 3.166112194e-04 0.021964041570 3.618904520e-08 7.4 309 3.452669830e-04 0.021727935579 3.049675551e-07 6.5 310 3.765163148e-04 0.021491829587 3.034669259e-07 6.5 311 4.105939528e-04 0.021255723595 6.315431433e-08 7.2 312 4.477558805e-04 0.021019617603 2.569936886e-08 7.6 313 4.882812500e-04 0.020783511612 -1.990711820e-07 6.7 314 5.324744788e-04 0.020547405620 -1.808471892e-07 6.7 315 5.806675366e-04 0.020311299628 9.470999129e-08 7.0 316 6.332224388e-04 0.020075193637 2.537703486e-07 6.6 317 6.905339660e-04 0.019839087645 -5.540736026e-08 7.3 318 7.530326296e-04 0.019602981653 -3.165362417e-08 7.5 319 8.211879055e-04 0.019366875662 -3.256026071e-08 7.5 320 8.955117609e-04 0.019130769670 -2.782856834e-08 7.6 321 9.765625000e-04 0.018894663678 1.266636318e-08 7.9 322 1.064948958e-03 0.018658557687 -2.358667284e-07 6.6 323 1.161335073e-03 0.018422451695 -2.867862725e-08 7.5 324 1.266444878e-03 0.018186345703 -4.015371813e-08 7.4 325 1.381067932e-03 0.017950239712 -2.243580433e-07 6.6 326 1.506065259e-03 0.017714133720 1.295179664e-07 6.9 327 1.642375811e-03 0.017478027728 5.258487723e-08 7.3 328 1.791023522e-03 0.017241921737 -7.191778106e-08 7.1 329 1.953125000e-03 0.017005815745 1.168092258e-07 6.9 330 2.129897915e-03 0.016769709753 2.579736935e-08 7.6 331 2.322670146e-03 0.016533603761 2.070165972e-08 7.7 332 2.532889755e-03 0.016297497770 1.453019420e-07 6.8 333 2.762135864e-03 0.016061391778 -1.448045621e-07 6.8 334 3.012130518e-03 0.015825285786 9.153762492e-08 7.0 335 3.284751622e-03 0.015589179795 3.166067630e-08 7.5 336 3.582047044e-03 0.015353073803 2.793761855e-08 7.6 337 3.906250000e-03 0.015116967811 -1.335779407e-07 6.9 338 4.259795831e-03 0.014880861820 -1.127023992e-07 6.9 339 4.645340293e-03 0.014644755828 1.472359789e-07 6.8 340 5.065779510e-03 0.014408649836 1.913901931e-07 6.7 341 5.524271728e-03 0.014172543845 1.078438131e-07 7.0 342 6.024261037e-03 0.013936437853 -1.618990075e-08 7.8 343 6.569503244e-03 0.013700331861 -9.444800253e-08 7.0 344 7.164094088e-03 0.013464225870 -4.169422585e-08 7.4 345 7.812500000e-03 0.013228119878 -1.890127832e-09 8.7 346 8.519591661e-03 0.012992013886 -9.842786630e-08 7.0 347 9.290680586e-03 0.012755907895 1.241692305e-10 9.9 348 1.013155902e-02 0.012519801903 -3.338975185e-08 7.5 349 1.104854346e-02 0.012283695911 1.347443089e-07 6.9 350 1.204852207e-02 0.012047589919 1.089463753e-08 8.0 351 1.313900649e-02 0.011811483928 -5.345811105e-08 7.3 352 1.432818818e-02 0.011575377936 7.911680988e-08 7.1 353 1.562500000e-02 0.011339271944 -1.082518675e-08 8.0 354 1.703918332e-02 0.011103165953 3.672092086e-08 7.4 355 1.858136117e-02 0.010867059961 3.491128742e-08 7.5 356 2.026311804e-02 0.010630953969 4.962150579e-09 8.3 357 2.209708691e-02 0.010394847978 -1.451724585e-08 7.8 358 2.409704415e-02 0.010158741986 3.193182674e-08 7.5 359 2.627801298e-02 0.009922635994 6.318548917e-08 7.2 360 2.865637635e-02 0.009686530003 3.477948718e-08 7.5 361 3.125000000e-02 0.009450424011 8.024744635e-08 7.1 362 3.407836665e-02 0.009214318019 7.269387325e-08 7.1 363 3.716272234e-02 0.008978212028 -6.037044309e-08 7.2 364 4.052623608e-02 0.008742106036 4.128745124e-08 7.4 365 4.419417382e-02 0.008506000044 -1.954348461e-09 8.7 366 4.819408829e-02 0.008269894053 6.703489763e-09 8.2 367 5.255602595e-02 0.008033788061 6.777183081e-08 7.2 368 5.731275270e-02 0.007797682069 3.394610371e-08 7.5 369 6.250000000e-02 0.007561576077 6.149548670e-08 7.2 370 6.815673329e-02 0.007325470086 -2.079005479e-08 7.7 371 7.432544469e-02 0.007089364094 2.085197826e-08 7.7 372 8.105247217e-02 0.006853258102 -3.090769152e-08 7.5 373 8.838834765e-02 0.006617152111 -5.460620156e-08 7.3 374 9.638817659e-02 0.006381046119 4.531310338e-08 7.3 375 1.051120519e-01 0.006144940127 1.717334919e-08 7.8 376 1.146255054e-01 0.005908834136 1.463540500e-08 7.8 377 1.250000000e-01 0.005672728144 3.113683844e-08 7.5 378 1.363134666e-01 0.005436622152 1.390943816e-08 7.9 379 1.486508894e-01 0.005200516161 -3.333184884e-08 7.5 380 1.621049443e-01 0.004964410169 -1.011861062e-08 8.0 381 1.767766953e-01 0.004728304177 -2.307868430e-08 7.6 382 1.927763532e-01 0.004492198186 -3.451609221e-09 8.5 383 2.102241038e-01 0.004256092194 -1.417570927e-08 7.8 384 2.292510108e-01 0.004019986202 -5.527056146e-09 8.3 385 2.500000000e-01 0.003783880211 -9.496292419e-09 8.0 386 2.726269332e-01 0.003547774219 -4.487596073e-09 8.3 387 2.973017788e-01 0.003311668227 -4.862096725e-09 8.3 388 3.242098887e-01 0.003075562235 -7.645436728e-09 8.1 389 3.535533906e-01 0.002839456244 3.773244051e-09 8.4 390 3.855527064e-01 0.002603350252 5.646113910e-09 8.2 391 4.204482076e-01 0.002367244260 -3.754701661e-09 8.4 392 4.585020216e-01 0.002131138269 -2.839019464e-09 8.5 393 5.000000000e-01 0.001895032277 -5.991975804e-10 9.2 > > ## } ## only when we find inaccurate regions > showProc.time() Time elapsed: 0.16 0.02 0.19 > > > ## 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: 8.279884e-15" > 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 elapsed: 0.03 0 0.06 > > 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 elapsed: 1.42 0.01 1.51 > > ## 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 elapsed: 4.09 0.05 4.17 > > ## 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 elapsed: 0.87 0.02 0.93 > > ##--- 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 elapsed: 0.25 0 0.25 > > 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 elapsed: 0 0 0 > > 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 elapsed: 0.77 0.03 0.86 > > ##-- 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 elapsed: 0.54 0 0.61 > > ## 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 elapsed: 0.11 0.01 0.13 > > ###------------- 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": > > if(!exists("lseq", mode="function")) + 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", lty=1, xlim=c(8e-9, 1e-3)) + abline(v= 3.47548562941137e-08, col = "gray80", lwd=3)#<- the cutoff value of lgamma1p() + + ## 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 32 bits per limb Attaching package: 'Rmpfr' The following object is masked from 'package:gmp': outer The following objects are masked from 'package:stats': dbinom, dgamma, dnorm, dpois, 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.1810533509579636008827292017839333436416956277420226420324997e-57" Time elapsed: 1.97 0 2.01 Class 'mpfrArray' [package "Rmpfr"] of dimension c(125, 7) and precision 200 -1 -1 ... num [1:125, 1:7] -1.00 -1.00 6.34e+13 3.64e+13 -1.00 ... - attr(*, "dimnames")=List of 2 ..$ : NULL ..$ : NULL Time elapsed: 0.08 0 0.09 Warning messages: 1: In xy.coords(x, y, xlabel, ylabel, log = log) : 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() approxmation > ## ----------------- --> .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 elapsed: 0.01 0 0.02 > > > ## 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 elapsed: 0.02 0 0.01 > 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: 3.820202e-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: 3.903776e-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] "Mean relative difference: 1.72789e-15" > ## 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) : 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) : 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) : 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) : 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) : 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) : 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) : 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) : 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) : 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) : 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) : 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) : 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) : 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) : 845 y values <= 0 omitted from logarithmic plot > summary(warnings()) 1 identical warnings: In xy.coords(x, y, xlabel, ylabel, log = log) : 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) : 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) : 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) : 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) : 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) : 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) : 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) : 226 y values <= 0 omitted from logarithmic plot > showProc.time() Time elapsed: 7.72 0.08 8.05 > > ## 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) : 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) : 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) : 363 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) : 363 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 elapsed: 0.11 0 0.11 > > ###--> 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 -1.27e-116 $ iter : int 8 $ init.it : int NA $ estim.prec: num 4.16e-17 > 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 elapsed: 0.11 0.02 0.16 > ## 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) : 125 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] 1.687539e-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) Warning message: In plot.window(...) : relative range of values = 90 * EPS, is small (axis 1) > > showProc.time() Time elapsed: 0.37 0.03 0.45 > > 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: 1, Adjusted R-squared: 1 F-statistic: 7.853e+13 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: 1.359e+16 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. 0.000e+00 6.328e-15 1.293e-14 1.575e-14 2.465e-14 3.741e-14 > plot(a, abs(1 - y/1e-100), log="xy", type="b") Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 1 y value <= 0 omitted from logarithmic plot > 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.302676e-16" > stopifnot(all.equal(z1[1:10], z1true, tol = 1e-15)) > showProc.time() Time elapsed: 0.17 0.03 0.2 > ##> > ##> 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") Error in axis(side = side, at = at, labels = labels, ...) : CreateAtVector [log-axis()]: axp[0] = 0 < 0! Calls: curve ... plot.default -> localAxis -> Axis -> Axis.default -> axis In addition: Warning messages: 1: In xy.coords(x, y, xlabel, ylabel, log) : 514 y values <= 0 omitted from logarithmic plot 2: In plot.window(...) : Internal(pretty()): very small range.. corrected 3: In axis(side = side, at = at, labels = labels, ...) : CreateAtVector "log"(from axis()): axp[0] = 0 ! Execution halted ** running tests for arch 'x64' ... [121s] ERROR Running 'chisq-nonc-ex.R' [33s] Running 'dnchisq-tst.R' [1s] Running 'pnbeta-tst.R' [1s] Running 'pnt-precision-problem.R' [21s] Running 'ppois-ex.R' [2s] Running 'qPoisBinom-ex.R' [1s] Running 'qbeta-dist.R' [15s] Running 'qbeta-tst.R' [1s] Running 'qgamma-ex.R' [24s] Running 't-nonc-tst.R' [10s] Running 'wienergerm-accuracy.R' [8s] Running 'wienergerm-pchisq-tst.R' [1s] Running 'wienergerm_nchisq.R' [3s] Running the tests in 'tests/qgamma-ex.R' failed. Complete output: > 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="Matrix", "test-tools-1.R", mustWork=TRUE)) Loading required package: tools > ##--> showProc.time(), assertError(), relErrV(), ... > showProc.time() Time elapsed: 1.85 0.06 1.98 > > (doExtras <- DPQ:::doExtras()) [1] FALSE > (sdir <- system.file("safe", package="DPQ")) ## save directory (to read from) [1] "D:/temp/RtmpcJZjAr/RLIBS_dec02f1a3c6f/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 elapsed: 0.83 0.05 0.89 > > require(graphics) > if(!dev.interactive(orNone=TRUE)) pdf("qgamma-appr.pdf") > eaxis <- sfsmisc::eaxis > > showProc.time() Time elapsed: 0.09 0.01 0.11 > ## 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 6.46593024678e-08 7.2 [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 3.94217658517e-10 9.4 [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 7.75564700239e-08 7.1 [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 -7.17404703909e-06 5.1 [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 -3.00271327148e-04 3.5 [20,] 0.014753716559535 0.01155624353203 1.52962087441e-10 9.8 [21,] 0.015499213434879 0.01125420134457 1.00492447767e-04 4.0 [22,] 0.015499213434879 0.01135920381800 -9.55739012376e-08 7.0 [23,] 0.018603016576955 0.01071716109330 -2.07537936112e-03 2.7 [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 4.43079705148e-08 7.4 [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 -1.78633419168e-09 8.7 [31,] 0.081818424963746 0.00689856335721 -2.30915286892e-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 -6.63889809793e-09 8.2 [36,] 0.110869751789691 0.00603268830251 8.95578944338e-10 9.0 [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 8.96455243371e-10 9.0 [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.78565421871e-06 5.7 [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.39812605937e-11 10.9 [59,] 0.214780951412088 0.00419869272965 -1.91879370171e-08 7.7 [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 2.14569900181e-07 6.7 [64,] 0.253082914021191 0.00375259362798 3.65436092498e-09 8.4 [65,] 0.253922058700076 0.00371237348323 -5.09014937222e-06 5.3 [66,] 0.254289278570932 0.00374343873151 -1.05664899053e-09 9.0 [67,] 0.260017499519858 0.00366179605930 -2.90430199223e-07 6.5 [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 -6.07109695849e-10 9.2 [72,] 0.284601233201101 0.00344936339590 -2.63026489478e-10 9.6 [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 -4.48927431229e-07 6.3 [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.82253101499e-11 10.7 [81,] 0.303182283998467 0.00328723648952 -2.23672191879e-11 10.7 [82,] 0.322319846303892 0.00306134512927 -1.15130830540e-05 4.9 [83,] 0.322319846303892 0.00310689001755 8.57751647487e-11 10.1 [84,] 0.325071272052651 0.00302343293053 -2.47088704493e-04 3.6 [85,] 0.325071272052651 0.00304146419577 -7.71376615361e-06 5.1 [86,] 0.331888412404218 0.00300837121343 -4.96098895297e-09 8.3 [87,] 0.362278153188527 0.00278204202032 4.53939330569e-10 9.3 [88,] 0.385389476781711 0.00260981704384 -1.77430203863e-09 8.8 [89,] 0.425333956955001 0.00232995789362 1.82823025607e-08 7.7 [90,] 0.439503709203564 0.00222452690840 -4.53585193982e-06 5.3 [91,] 0.439503709203564 0.00224964327069 -3.02331937263e-10 9.5 [92,] 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 6.465930e-08 7.2 [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 7.755647e-08 7.1 [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 -7.174047e-06 5.1 [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 -3.002713e-04 3.5 [17,] 0.0154992134 0.011254201 1.004924e-04 4.0 [18,] 0.0154992134 0.011359204 -9.557390e-08 7.0 [19,] 0.0186030166 0.010717161 -2.075379e-03 2.7 [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 4.430797e-08 7.4 [23,] 0.0237302174 0.010162521 -1.077327e-06 6.0 [24,] 0.0447535254 0.008396264 1.222242e-05 4.9 [25,] 0.0818184250 0.006860077 -1.786334e-09 8.7 [26,] 0.0828003091 0.006812347 4.179976e-09 8.4 [27,] 0.0908216581 0.006552698 -7.160336e-09 8.1 [28,] 0.1022947605 0.006235631 -6.638898e-09 8.2 [29,] 0.1274058577 0.005623691 6.605415e-09 8.2 [30,] 0.1352296342 0.005400734 -2.347626e-05 4.6 [31,] 0.1377322800 0.005330921 2.992858e-04 3.5 [32,] 0.1381129175 0.005351387 -2.053358e-06 5.7 [33,] 0.1411006360 0.005273058 4.315938e-05 4.4 [34,] 0.1429052994 0.005236800 3.481808e-04 3.5 [35,] 0.1456245572 0.005269240 -1.945018e-09 8.7 [36,] 0.1546068729 0.005068069 -4.599247e-07 6.3 [37,] 0.1546068729 0.005073662 2.723010e-07 6.6 [38,] 0.1697410365 0.004841818 5.314010e-09 8.3 [39,] 0.1773275763 0.004659561 5.534044e-05 4.3 [40,] 0.1900940174 0.004503736 1.785654e-06 5.7 [41,] 0.1901476415 0.004534687 5.662356e-09 8.2 [42,] 0.2147809514 0.004198693 -1.918794e-08 7.7 [43,] 0.2328051066 0.003953993 -9.175810e-08 7.0 [44,] 0.2520765119 0.003749038 2.145699e-07 6.7 [45,] 0.2530829140 0.003752594 3.654361e-09 8.4 [46,] 0.2539220587 0.003712373 -5.090149e-06 5.3 [47,] 0.2542892786 0.003743439 -1.056649e-09 9.0 [48,] 0.2600174995 0.003661796 -2.904302e-07 6.5 [49,] 0.2703239068 0.003519992 -1.561648e-04 3.8 [50,] 0.2716993561 0.003550681 5.130930e-09 8.3 [51,] 0.2755161961 0.003468040 -4.351715e-04 3.4 [52,] 0.2901885431 0.003366135 -5.644431e-08 7.2 [53,] 0.2905790220 0.003344235 1.026676e-07 7.0 [54,] 0.2905790220 0.003367649 2.260616e-08 7.6 [55,] 0.2918501987 0.003335528 -1.273388e-06 5.9 [56,] 0.2965211365 0.003303089 -4.489274e-07 6.3 [57,] 0.2980341749 0.003304623 8.424704e-09 8.1 [58,] 0.3005567833 0.003239225 4.660033e-05 4.3 [59,] 0.3223198463 0.003061345 -1.151308e-05 4.9 [60,] 0.3250712721 0.003023433 -2.470887e-04 3.6 [61,] 0.3250712721 0.003041464 -7.713766e-06 5.1 [62,] 0.3318884124 0.003008371 -4.960989e-09 8.3 [63,] 0.3853894768 0.002609817 -1.774302e-09 8.8 [64,] 0.4253339570 0.002329958 1.828230e-08 7.7 [65,] 0.4395037092 0.002224527 -4.535852e-06 5.3 [66,] 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.912e-04 -1.441e-04 -1.969e-05 7.710e-05 1.082e-03 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 6.184e-06 1.133e-05 0.546 0.587 log(p) -2.725e-03 3.660e-06 -744.530 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.0002557 on 90 degrees of freedom Multiple R-squared: 0.9998, Adjusted R-squared: 0.9998 F-statistic: 5.543e+05 on 1 and 90 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.094448579669 -5.693322411e-07 6.2 2 9.685654347e-16 0.094212473681 7.322734105e-07 6.1 3 1.056228096e-15 0.093976367693 -5.153676885e-08 7.3 4 1.151824906e-15 0.093740261705 -8.032222418e-07 6.1 5 1.256073967e-15 0.093504155717 5.733035913e-07 6.2 6 1.369758374e-15 0.093268049728 -1.195090105e-07 6.9 7 1.493732098e-15 0.093031943740 -7.801879547e-07 6.1 8 1.628926404e-15 0.092795837752 6.712882394e-07 6.2 9 1.776356839e-15 0.092559731764 6.950104159e-08 7.2 10 1.937130869e-15 0.092323625776 -5.001435726e-07 6.3 11 2.112456192e-15 0.092087519788 1.026313132e-06 6.0 12 2.303649813e-15 0.091851413800 -1.542993703e-06 5.8 13 2.512147934e-15 0.091615307811 3.699662154e-08 7.4 14 2.739516748e-15 0.091379201823 -4.094315547e-07 6.4 15 2.987464197e-15 0.091143095835 -8.237025693e-07 6.1 16 3.257852808e-15 0.090906989847 8.313182130e-07 6.1 17 3.552713679e-15 0.090670883859 4.759887695e-07 6.3 18 3.874261739e-15 0.090434777871 1.528253788e-07 6.8 19 4.224912384e-15 0.090198671883 -1.381691319e-07 6.9 20 4.607299625e-15 0.089962565895 -3.969919471e-07 6.4 21 5.024295868e-15 0.089726459906 1.386675560e-06 5.9 22 5.479033495e-15 0.089490353918 -8.181112241e-07 6.1 23 5.974928394e-15 0.089254247930 1.019155402e-06 6.0 24 6.515705616e-15 0.089018141942 -1.110509928e-06 6.0 25 7.105427358e-15 0.088782035954 7.803687333e-07 6.1 26 7.748523477e-15 0.088545929966 -1.274165560e-06 5.9 27 8.449824769e-15 0.088309823978 6.703380386e-07 6.2 28 9.214599250e-15 0.088073717989 -1.309055653e-06 5.9 29 1.004859174e-14 0.087837612001 6.890857738e-07 6.2 30 1.095806699e-14 0.087601506013 -1.215157786e-06 5.9 31 1.194985679e-14 0.087365400025 8.366343469e-07 6.1 32 1.303141123e-14 0.087129294037 -9.924495552e-07 6.0 33 1.421085472e-14 0.086893188049 1.113006151e-06 6.0 34 1.549704695e-14 0.086657082061 -6.409085871e-07 6.2 35 1.689964954e-14 0.086420976073 -4.168188492e-07 6.4 36 1.842919850e-14 0.086184870084 -1.605125546e-07 6.8 37 2.009718347e-14 0.085948764096 1.280130847e-07 6.9 38 2.191613398e-14 0.085712658108 4.487608530e-07 6.3 39 2.389971358e-14 0.085476552120 -1.111735280e-06 6.0 40 2.606282246e-14 0.085240446132 1.186933915e-06 5.9 41 2.842170943e-14 0.085004340144 -2.983612493e-07 6.5 42 3.099409391e-14 0.084768234156 1.566736150e-07 6.8 43 3.379929908e-14 0.084532128167 6.439440862e-07 6.2 44 3.685839700e-14 0.084296022179 -7.230811827e-07 6.1 45 4.019436694e-14 0.084059916191 -1.659619306e-07 6.8 46 4.383226796e-14 0.083823810203 4.234008103e-07 6.4 47 4.779942715e-14 0.083587704215 -8.253375789e-07 6.1 48 5.212564492e-14 0.083351598227 -1.661117941e-07 6.8 49 5.684341886e-14 0.083115492239 5.253653413e-07 6.3 50 6.198818782e-14 0.082879386251 -6.050692618e-07 6.2 51 6.759859815e-14 0.082643280262 1.562851194e-07 6.8 52 7.371679400e-14 0.082407174274 9.498987034e-07 6.0 53 8.038873388e-14 0.082171068286 -6.221528537e-08 7.2 54 8.766453592e-14 0.081934962298 -1.031256456e-06 6.0 55 9.559885430e-14 0.081698856310 -1.301203456e-07 6.9 56 1.042512898e-13 0.081462750322 8.032851107e-07 6.1 57 1.136868377e-13 0.081226644334 -4.741459603e-08 7.3 58 1.239763756e-13 0.080990538345 -8.550366157e-07 6.1 59 1.351971963e-13 0.080754432357 1.859174034e-07 6.7 60 1.474335880e-13 0.080518326369 1.259150868e-06 5.9 61 1.607774678e-13 0.080282220381 5.698910264e-07 6.2 62 1.753290718e-13 0.080046114393 -7.628630283e-08 7.1 63 1.911977086e-13 0.079810008405 -6.793803771e-07 6.2 64 2.085025797e-13 0.079573902417 -1.239390442e-06 5.9 65 2.273736754e-13 0.079337796429 1.671630268e-08 7.8 66 2.479527513e-13 0.079101690440 1.305117147e-06 5.9 67 2.703943926e-13 0.078865584452 8.634978415e-07 6.1 68 2.948671760e-13 0.078629478464 4.649672880e-07 6.3 69 3.215549355e-13 0.078393372476 1.095262072e-07 7.0 70 3.506581437e-13 0.078157266488 -2.028246688e-07 6.7 71 3.823954172e-13 0.077921160500 -4.720846172e-07 6.3 72 4.170051594e-13 0.077685054512 -6.982529182e-07 6.2 73 4.547473509e-13 0.077448948523 -8.813288566e-07 6.1 74 4.959055026e-13 0.077212842535 -1.021311716e-06 6.0 75 5.407887852e-13 0.076976736547 -1.118200798e-06 6.0 76 5.897343520e-13 0.076740630559 -1.171995389e-06 5.9 77 6.431098711e-13 0.076504524571 -1.182694798e-06 5.9 78 7.013162874e-13 0.076268418583 -1.150298324e-06 5.9 79 7.647908344e-13 0.076032312595 -1.074805275e-06 6.0 80 8.340103188e-13 0.075796206607 -9.562149603e-07 6.0 81 9.094947018e-13 0.075560100618 -7.945266962e-07 6.1 82 9.918110051e-13 0.075323994630 -5.897397957e-07 6.2 83 1.081577570e-12 0.075087888642 -3.418535914e-07 6.5 84 1.179468704e-12 0.074851782654 -5.086740451e-08 7.3 85 1.286219742e-12 0.074615676666 2.832194348e-07 6.5 86 1.402632575e-12 0.074379570678 6.604075928e-07 6.2 87 1.529581669e-12 0.074143464690 1.080697732e-06 6.0 88 1.668020638e-12 0.073907358701 -1.042647599e-07 7.0 89 1.818989404e-12 0.073671252713 4.076445904e-07 6.4 90 1.983622010e-12 0.073435146725 -6.748022854e-07 6.2 91 2.163155141e-12 0.073199040737 -7.127495261e-08 7.1 92 2.358937408e-12 0.072962934749 5.753558137e-07 6.2 93 2.572439484e-12 0.072726828761 -3.560679425e-07 6.4 94 2.805265149e-12 0.072490722773 3.821799130e-07 6.4 95 3.059163338e-12 0.072254616785 -4.467414887e-07 6.3 96 3.336041275e-12 0.072018510796 3.831221452e-07 6.4 97 3.637978807e-12 0.071782404808 -3.433027738e-07 6.5 98 3.967244020e-12 0.071546298820 -1.015745335e-06 6.0 99 4.326310282e-12 0.071310192832 -4.575907031e-08 7.3 100 4.717874816e-12 0.071074086844 9.673320129e-07 6.0 101 5.144878969e-12 0.070837980856 -1.131696160e-06 5.9 102 5.610530299e-12 0.070601874868 -2.159353318e-08 7.7 103 6.118326675e-12 0.070365768879 1.131613904e-06 5.9 104 6.672082550e-12 0.070129662891 -7.946347651e-07 6.1 105 7.275957614e-12 0.069893556903 4.555778093e-07 6.3 106 7.934488041e-12 0.069657450915 1.985071415e-07 6.7 107 8.652620563e-12 0.069421344927 -4.602383807e-09 8.3 108 9.435749632e-12 0.069185238939 -1.537536061e-07 6.8 109 1.028975794e-11 0.068949132951 -2.489493711e-07 6.6 110 1.122106060e-11 0.068713026963 -2.901925258e-07 6.5 111 1.223665335e-11 0.068476920974 -2.774859360e-07 6.6 112 1.334416510e-11 0.068240814986 -2.108324577e-07 6.7 113 1.455191523e-11 0.068004708998 -9.023495839e-08 7.0 114 1.586897608e-11 0.067768603010 8.430368981e-08 7.1 115 1.730524113e-11 0.067532497022 3.127806092e-07 6.5 116 1.887149926e-11 0.067296391034 -9.005818178e-07 6.0 117 2.057951587e-11 0.067060285046 9.315377304e-07 6.0 118 2.244212120e-11 0.066824179057 -1.630697595e-07 6.8 119 2.447330670e-11 0.066588073069 2.865757222e-07 6.5 120 2.668833020e-11 0.066351967081 7.901436261e-07 6.1 121 2.910383046e-11 0.066115861093 -1.208588474e-07 6.9 122 3.173795216e-11 0.065879755105 -9.670019019e-07 6.0 123 3.461048225e-11 0.065643649117 -2.908072072e-07 6.5 124 3.774299853e-11 0.065407543129 4.392954942e-07 6.4 125 4.115903175e-11 0.065171437141 -2.233106784e-07 6.7 126 4.488424239e-11 0.064935331152 -8.210848088e-07 6.1 127 4.894661340e-11 0.064699225164 8.158930909e-08 7.1 128 5.337666040e-11 0.064463119176 1.038156913e-06 6.0 129 5.820766091e-11 0.064227013188 6.238512490e-07 6.2 130 6.347590433e-11 0.063990907200 2.743501614e-07 6.6 131 6.922096451e-11 0.063754801212 -1.035416930e-08 8.0 132 7.548599706e-11 0.063518695224 -2.302695521e-07 6.6 133 8.231806350e-11 0.063282589235 -3.854038086e-07 6.4 134 8.976848478e-11 0.063046483247 -4.757647607e-07 6.3 135 9.789322680e-11 0.062810377259 -5.013602509e-07 6.3 136 1.067533208e-10 0.062574271271 -4.621981136e-07 6.3 137 1.164153218e-10 0.062338165283 -3.582861912e-07 6.4 138 1.269518087e-10 0.062102059295 -1.896323394e-07 6.7 139 1.384419290e-10 0.061865953307 4.375558904e-08 7.4 140 1.509719941e-10 0.061629847319 3.418697269e-07 6.5 141 1.646361270e-10 0.061393741330 7.047022075e-07 6.2 142 1.795369696e-10 0.061157635342 -2.212484238e-07 6.7 143 1.957864536e-10 0.060921529354 2.764617620e-07 6.6 144 2.135066416e-10 0.060685423366 -5.036437378e-07 6.3 145 2.328306437e-10 0.060449317378 1.289049596e-07 6.9 146 2.539036173e-10 0.060213211390 8.261288568e-07 6.1 147 2.768838580e-10 0.059977105402 2.619449679e-07 6.6 148 3.019439882e-10 0.059740999413 -2.266130656e-07 6.6 149 3.292722540e-10 0.059504893425 6.754948677e-07 6.2 150 3.590739391e-10 0.059268787437 -9.769083436e-07 6.0 151 3.915729072e-10 0.059032681449 6.536886865e-08 7.2 152 4.270132832e-10 0.058796575461 -1.263299489e-07 6.9 153 4.656612873e-10 0.058560469473 -2.424715895e-07 6.6 154 5.078072346e-10 0.058324363485 -2.830702768e-07 6.5 155 5.537677160e-10 0.058088257497 -2.481402319e-07 6.6 156 6.038879765e-10 0.057852151508 -1.376956813e-07 6.9 157 6.585445080e-10 0.057616045520 4.824913535e-08 7.3 158 7.181478783e-10 0.057379939532 3.096799803e-07 6.5 159 7.831458144e-10 0.057143833544 6.465826115e-07 6.2 160 8.540265665e-10 0.056907727556 -1.956638767e-07 6.7 161 9.313225746e-10 0.056671621568 2.976208309e-07 6.5 162 1.015614469e-09 0.056435515580 8.663322291e-07 6.1 163 1.107535432e-09 0.056199409591 2.723392490e-07 6.6 164 1.207775953e-09 0.055963303603 -2.352547019e-07 6.6 165 1.317089016e-09 0.055727197615 -6.564715971e-07 6.2 166 1.436295757e-09 0.055491091627 2.302113088e-07 6.6 167 1.566291629e-09 0.055254985639 -2.383583730e-08 7.6 168 1.708053133e-09 0.055018879651 -1.915690806e-07 6.7 169 1.862645149e-09 0.054782773663 -2.730104096e-07 6.6 170 2.031228938e-09 0.054546667675 -2.681818159e-07 6.6 171 2.215070864e-09 0.054310561686 -1.771053018e-07 6.8 172 2.415551906e-09 0.054074455698 1.971357522e-10 9.7 173 2.634178032e-09 0.053838349710 2.637034958e-07 6.6 174 2.872591513e-09 0.053602243722 -5.640842873e-07 6.2 175 3.132583258e-09 0.053366137734 -1.227401389e-07 6.9 176 3.416106266e-09 0.053130031746 4.047391041e-07 6.4 177 3.725290298e-09 0.052893925758 -1.426154970e-07 6.8 178 4.062457877e-09 0.052657819769 -5.928550000e-07 6.2 179 4.430141728e-09 0.052421713781 2.038662528e-07 6.7 180 4.831103812e-09 0.052185607793 -5.776660905e-08 7.2 181 5.268356064e-09 0.051949501805 -2.223745084e-07 6.7 182 5.745183026e-09 0.051713395817 8.433560312e-07 6.1 183 6.265166516e-09 0.051477289829 -2.606399732e-07 6.6 184 6.832212532e-09 0.051241183841 -1.343598028e-07 6.9 185 7.450580597e-09 0.051005077853 8.882078772e-08 7.1 186 8.124915754e-09 0.050768971864 -7.023640207e-07 6.2 187 8.860283457e-09 0.050532865876 8.257587175e-07 6.1 188 9.662207623e-09 0.050296759888 2.392302341e-07 6.6 189 1.053671213e-08 0.050060653900 -2.394743450e-07 6.6 190 1.149036605e-08 0.049824547912 -6.103966785e-07 6.2 191 1.253033303e-08 0.049588441924 2.100143996e-07 6.7 192 1.366442506e-08 0.049352335936 4.899553085e-08 7.3 193 1.490116119e-08 0.049116229947 -4.362272321e-09 8.4 194 1.624983151e-08 0.048880123959 4.989935698e-08 7.3 195 1.772056691e-08 0.048644017971 2.117388018e-07 6.7 196 1.932441525e-08 0.048407911983 -5.748355087e-07 6.2 197 2.107342426e-08 0.048171805995 -1.924479196e-07 6.7 198 2.298073210e-08 0.047935700007 2.973889422e-07 6.5 199 2.506066606e-08 0.047699594019 -1.447314699e-07 6.8 200 2.732885013e-08 0.047463488031 -4.684300434e-07 6.3 201 2.980232239e-08 0.047227382042 3.545092193e-07 6.5 202 3.249966302e-08 0.046991276054 -7.607756496e-07 6.1 203 3.544113383e-08 0.046755170066 2.876612886e-07 6.5 204 3.864883049e-08 0.046519064078 -5.800756444e-07 6.2 205 4.214684851e-08 0.046282958090 6.936620679e-07 6.2 206 4.596146421e-08 0.046046852102 7.324251272e-08 7.1 207 5.012133212e-08 0.045810746114 -4.180421775e-07 6.4 208 5.465770025e-08 0.045574640125 2.092254814e-07 6.7 209 5.960464478e-08 0.045338534137 -2.954279021e-08 7.5 210 6.499932603e-08 0.045102428149 -1.393715623e-07 6.9 211 7.088226765e-08 0.044866322161 -1.203274604e-07 6.9 212 7.729766099e-08 0.044630216173 2.752292738e-08 7.6 213 8.429369702e-08 0.044394110185 -6.576521183e-07 6.2 214 9.192292842e-08 0.044158004197 -2.468590250e-07 6.6 215 1.002426642e-07 0.043921898209 2.925361181e-07 6.5 216 1.093154005e-07 0.043685792220 1.531700988e-08 7.8 217 1.192092896e-07 0.043449686232 -1.223621087e-07 6.9 218 1.299986521e-07 0.043213580244 -1.205823330e-07 6.9 219 1.417645353e-07 0.042977474256 2.057528392e-08 7.7 220 1.545953220e-07 0.042741368268 -6.218908339e-07 6.2 221 1.685873940e-07 0.042505262280 -1.966860006e-07 6.7 222 1.838458568e-07 0.042269156292 3.676487106e-07 6.4 223 2.004853285e-07 0.042033050303 1.647381379e-07 6.8 224 2.186308010e-07 0.041796944315 1.118715082e-07 7.0 225 2.384185791e-07 0.041560838327 2.089520004e-07 6.7 226 2.599973041e-07 0.041324732339 4.558828597e-07 6.3 227 2.835290706e-07 0.041088626351 -3.149587724e-08 7.5 228 3.091906439e-07 0.040852520363 -3.581121097e-07 6.4 229 3.371747881e-07 0.040616414375 3.488470744e-07 6.5 230 3.676917137e-07 0.040380308387 -5.295130916e-07 6.3 231 4.009706570e-07 0.040144202398 4.872893066e-07 6.3 232 4.372616020e-07 0.039908096410 -5.923112045e-08 7.2 233 4.768371582e-07 0.039671990422 -4.344346312e-07 6.4 234 5.199946082e-07 0.039435884434 2.066676108e-07 6.7 235 5.670581412e-07 0.039199778446 1.681372633e-07 6.8 236 6.183812879e-07 0.038963672458 -5.334671598e-07 6.3 237 6.743495762e-07 0.038727566470 -2.247258435e-07 6.6 238 7.353834273e-07 0.038491460481 2.546716327e-07 6.6 239 8.019413140e-07 0.038255354493 8.725913336e-08 7.1 240 8.745232040e-07 0.038019248505 1.013349561e-07 7.0 241 9.536743164e-07 0.037783142517 -5.094306235e-07 6.3 242 1.039989216e-06 0.037547036529 -1.272824359e-07 6.9 243 1.134116282e-06 0.037310930541 4.358925881e-07 6.4 244 1.236762576e-06 0.037074824553 3.904291315e-07 6.4 245 1.348699152e-06 0.036838718565 5.367821070e-07 6.3 246 1.470766855e-06 0.036602612576 9.641457999e-08 7.0 247 1.603882628e-06 0.036366506588 -1.413544166e-07 6.8 248 1.749046408e-06 0.036130400600 -1.767197404e-07 6.8 249 1.907348633e-06 0.035894294612 -9.875985141e-09 8.0 250 2.079978433e-06 0.035658188624 -3.970802771e-07 6.4 251 2.268232565e-06 0.035422082636 1.791465434e-07 6.7 252 2.473525152e-06 0.035185976648 -5.328653061e-07 6.3 253 2.697398305e-06 0.034949870659 -2.818818945e-07 6.5 254 2.941533709e-06 0.034713764671 1.813954295e-07 6.7 255 3.207765256e-06 0.034477658683 1.285281512e-07 6.9 256 3.498092816e-06 0.034241552695 2.986249925e-07 6.5 257 3.814697266e-06 0.034005446707 -2.562239221e-08 7.6 258 4.159956866e-06 0.033769340719 -1.162612726e-07 6.9 259 4.536465130e-06 0.033533234731 2.644014629e-08 7.6 260 4.947050303e-06 0.033297128743 4.022141099e-07 6.4 261 5.394796609e-06 0.033061022754 -3.787335339e-07 6.4 262 5.883067419e-06 0.032824916766 4.734672679e-07 6.3 263 6.415530512e-06 0.032588810778 1.906558875e-07 6.7 264 6.996185632e-06 0.032352704790 1.620354394e-07 6.8 265 7.629394531e-06 0.032116598802 3.872828165e-07 6.4 266 8.319913732e-06 0.031880492814 -4.676686218e-07 6.3 267 9.072930260e-06 0.031644386826 2.754435448e-07 6.6 268 9.894100606e-06 0.031408280837 -4.002658072e-08 7.4 269 1.078959322e-05 0.031172174849 -8.071856361e-08 7.1 270 1.176613484e-05 0.030936068861 1.529846981e-07 6.8 271 1.283106102e-05 0.030699962873 2.171783053e-08 7.7 272 1.399237126e-05 0.030463856885 1.752040166e-07 6.8 273 1.525878906e-05 0.030227750897 -1.479256340e-08 7.8 274 1.663982746e-05 0.029991644909 9.026352044e-08 7.0 275 1.814586052e-05 0.029755538921 -1.267311944e-07 6.9 276 1.978820121e-05 0.029519432932 -3.843675112e-08 7.4 277 2.157918644e-05 0.029283326944 -2.508241883e-07 6.6 278 2.353226967e-05 0.029047220956 -1.477539424e-07 6.8 279 2.566212205e-05 0.028811114968 2.702582267e-07 6.6 280 2.798474253e-05 0.028575008980 -1.748270420e-07 6.8 281 3.051757812e-05 0.028338902992 -2.837500699e-07 6.5 282 3.327965493e-05 0.028102797004 -5.710114848e-08 7.2 283 3.629172104e-05 0.027866691015 -6.748044701e-08 7.2 284 3.957640242e-05 0.027630585027 2.676487311e-07 6.6 285 4.315837288e-05 0.027394479039 3.867922076e-07 6.4 286 4.706453935e-05 0.027158373051 -2.492501328e-07 6.6 287 5.132424410e-05 0.026922267063 4.134786513e-08 7.4 288 5.596948506e-05 0.026686161075 -3.911749389e-07 6.4 289 6.103515625e-05 0.026450055087 1.014410593e-07 7.0 290 6.655930986e-05 0.026213949098 -9.713157123e-08 7.0 291 7.258344208e-05 0.025977843110 1.004778277e-07 7.0 292 7.915280485e-05 0.025741737122 -3.501362080e-07 6.5 293 8.631674575e-05 0.025505631134 -3.839952303e-07 6.4 294 9.412907870e-05 0.025269525146 -2.027959578e-09 8.7 295 1.026484882e-04 0.025033419158 -2.152639500e-07 6.7 296 1.119389701e-04 0.024797313170 7.678417235e-09 8.1 297 1.220703125e-04 0.024561207182 -3.220049949e-07 6.5 298 1.331186197e-04 0.024325101193 -1.953370445e-07 6.7 299 1.451668842e-04 0.024088995205 -9.617589081e-08 7.0 300 1.583056097e-04 0.023852889217 -8.977143029e-09 8.0 301 1.726334915e-04 0.023616783229 8.175071808e-08 7.1 302 1.882581574e-04 0.023380677241 1.914434747e-07 6.7 303 2.052969764e-04 0.023144571253 -1.249542994e-07 6.9 304 2.238779402e-04 0.022908465265 7.430447435e-08 7.1 305 2.441406250e-04 0.022672359276 -1.108169194e-07 7.0 306 2.662372394e-04 0.022436253288 2.389977569e-07 6.6 307 2.903337683e-04 0.022200147300 2.458562984e-07 6.6 308 3.166112194e-04 0.021964041312 -5.850366769e-08 7.2 309 3.452669830e-04 0.021727935324 2.115776949e-07 6.7 310 3.765163148e-04 0.021491829336 2.113798955e-07 6.7 311 4.105939528e-04 0.021255723348 -2.762990770e-08 7.6 312 4.477558805e-04 0.021019617360 -6.378203188e-08 7.2 313 4.882812500e-04 0.020783511371 -2.872497766e-07 6.5 314 5.324744788e-04 0.020547405383 -2.677229627e-07 6.6 315 5.806675366e-04 0.020311299395 9.137057333e-09 8.0 316 6.332224388e-04 0.020075193407 1.695002406e-07 6.8 317 6.905339660e-04 0.019839087419 -1.383746775e-07 6.9 318 7.530326296e-04 0.019602981431 2.636260105e-07 6.6 319 8.211879055e-04 0.019366875443 -1.129219587e-07 6.9 320 8.955117609e-04 0.019130769454 2.589505672e-07 6.6 321 9.765625000e-04 0.018894663466 -6.508972539e-08 7.2 322 1.064948958e-03 0.018658557478 4.241218476e-08 7.4 323 1.161335073e-03 0.018422451490 2.453652983e-07 6.6 324 1.266444878e-03 0.018186345502 2.296493848e-07 6.6 325 1.381067932e-03 0.017950239514 4.120031982e-08 7.4 326 1.506065259e-03 0.017714133526 5.827584493e-08 7.2 327 1.642375811e-03 0.017478027538 -1.735446653e-08 7.8 328 1.791023522e-03 0.017241921549 1.809493316e-07 6.7 329 1.953125000e-03 0.017005815561 4.947543719e-08 7.3 330 2.129897915e-03 0.016769709573 -4.023365396e-08 7.4 331 2.322670146e-03 0.016533603585 -4.402659903e-08 7.4 332 2.532889755e-03 0.016297497597 8.187645184e-08 7.1 333 2.762135864e-03 0.016061391609 -2.069273115e-07 6.7 334 3.012130518e-03 0.015825285621 3.071764043e-08 7.5 335 3.284751622e-03 0.015589179632 -2.785656417e-08 7.6 336 3.582047044e-03 0.015353073644 -3.027688034e-08 7.5 337 3.906250000e-03 0.015116967656 -1.904897100e-07 6.7 338 4.259795831e-03 0.014880861668 -1.683114372e-07 6.8 339 4.645340293e-03 0.014644755680 9.292967829e-08 7.0 340 5.065779510e-03 0.014408649692 1.383866128e-07 6.9 341 5.524271728e-03 0.014172543704 5.614294096e-08 7.3 342 6.024261037e-03 0.013936437716 -6.658807195e-08 7.2 343 6.569503244e-03 0.013700331727 -1.435434784e-07 6.8 344 7.164094088e-03 0.013464225739 -8.948700758e-08 7.0 345 7.812500000e-03 0.013228119751 -4.838022560e-08 7.3 346 8.519591661e-03 0.012992013763 -1.436152919e-07 6.8 347 9.290680586e-03 0.012755907775 -4.376058627e-08 7.4 348 1.013155902e-02 0.012519801787 1.358225326e-07 6.9 349 1.104854346e-02 0.012283695799 9.346486218e-08 7.0 350 1.204852207e-02 0.012047589810 -2.908217578e-08 7.5 351 1.313900649e-02 0.011811483822 -9.213229890e-08 7.0 352 1.432818818e-02 0.011575377834 4.174524038e-08 7.4 353 1.562500000e-02 0.011339271846 1.379139563e-07 6.9 354 1.703918332e-02 0.011103165858 1.954543860e-09 8.7 355 1.858136117e-02 0.010867059870 1.447486908e-09 8.8 356 2.026311804e-02 0.010630953882 -2.719908920e-08 7.6 357 2.209708691e-02 0.010394847894 1.180522117e-07 6.9 358 2.409704415e-02 0.010158741905 2.375664332e-09 8.6 359 2.627801298e-02 0.009922635917 3.493183931e-08 7.5 360 2.865637635e-02 0.009686529929 7.828329762e-09 8.1 361 3.125000000e-02 0.009450423941 5.459876207e-08 7.3 362 3.407836665e-02 0.009214317953 4.834763923e-08 7.3 363 3.716272234e-02 0.008978211965 -8.341425239e-08 7.1 364 4.052623608e-02 0.008742105977 1.954604523e-08 7.7 365 4.419417382e-02 0.008505999988 -2.239338204e-08 7.6 366 4.819408829e-02 0.008269894000 -1.243319980e-08 7.9 367 5.255602595e-02 0.008033788012 4.993745217e-08 7.3 368 5.731275270e-02 0.007797682024 1.741400035e-08 7.8 369 6.250000000e-02 0.007561576036 4.626561989e-08 7.3 370 6.815673329e-02 0.007325470048 -3.471772847e-08 7.5 371 7.432544469e-02 0.007089364060 8.226451520e-09 8.1 372 8.105247217e-02 0.006853258072 -4.223112349e-08 7.4 373 8.838834765e-02 0.006617152083 1.657740989e-08 7.8 374 9.638817659e-02 0.006381046095 3.659368508e-08 7.4 375 1.051120519e-01 0.006144940107 9.755834918e-09 8.0 376 1.146255054e-01 0.005908834119 8.519715378e-09 8.1 377 1.250000000e-01 0.005672728131 2.632288454e-08 7.6 378 1.363134666e-01 0.005436622143 1.039711917e-08 8.0 379 1.486508894e-01 0.005200516155 -3.554264660e-08 7.4 380 1.621049443e-01 0.004964410166 -1.102801761e-08 8.0 381 1.767766953e-01 0.004728304178 -2.268685084e-08 7.6 382 1.927763532e-01 0.004492198190 -1.758708290e-09 8.8 383 2.102241038e-01 0.004256092202 -1.118194382e-08 8.0 384 2.292510108e-01 0.004019986214 -1.232663527e-09 8.9 385 2.500000000e-01 0.003783880226 -3.901555301e-09 8.4 386 2.726269332e-01 0.003547774238 2.407146371e-09 8.6 387 2.973017788e-01 0.003311668250 3.332240928e-09 8.5 388 3.242098887e-01 0.003075562261 1.847989872e-09 8.7 389 3.535533906e-01 0.002839456273 5.534314118e-10 9.3 390 3.855527064e-01 0.002603350285 -4.492217709e-09 8.3 391 4.204482076e-01 0.002367244297 1.137717809e-09 8.9 392 4.585020216e-01 0.002131138309 -5.251989954e-10 9.3 393 5.000000000e-01 0.001895032321 -1.494597335e-09 8.8 > > ## } ## only when we find inaccurate regions > showProc.time() Time elapsed: 0.13 0 0.15 > > > ## 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 elapsed: 0.03 0 0.06 > > 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 elapsed: 1.5 0 1.5 > > ## 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 elapsed: 4.41 0.06 4.57 > > ## 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 elapsed: 0.77 0.04 0.83 > > ##--- 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 elapsed: 0.15 0.01 0.19 > > 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 elapsed: 0 0 0 > > 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 elapsed: 0.6 0.05 0.7 > > ##-- 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 elapsed: 0.53 0 0.59 > > ## 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 elapsed: 0.07 0 0.11 > > ###------------- 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": > > if(!exists("lseq", mode="function")) + 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", lty=1, xlim=c(8e-9, 1e-3)) + abline(v= 3.47548562941137e-08, col = "gray80", lwd=3)#<- the cutoff value of lgamma1p() + + ## 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 objects are masked from 'package:stats': dbinom, dgamma, dnorm, dpois, 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.1810533509579636008827292017839333436416956277420226420324997e-57" Time elapsed: 1.39 0.05 1.48 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 elapsed: 0.1 0 0.1 Warning messages: 1: In xy.coords(x, y, xlabel, ylabel, log = log) : 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() approxmation > ## ----------------- --> .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 elapsed: 0 0 0.01 > > > ## 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 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) : 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) : 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) : 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) : 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) : 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) : 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) : 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) : 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) : 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) : 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) : 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) : 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) : 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) : 845 y values <= 0 omitted from logarithmic plot > summary(warnings()) 1 identical warnings: In xy.coords(x, y, xlabel, ylabel, log = log) : 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) : 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) : 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) : 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) : 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) : 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) : 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) : 226 y values <= 0 omitted from logarithmic plot > showProc.time() Time elapsed: 8.93 0.03 9.13 > > ## 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) : 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) : 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) : 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) : 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 elapsed: 0.07 0.01 0.09 > > ###--> 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 elapsed: 0.11 0 0.16 > ## 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) : 6 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) Warning message: In plot.window(...) : relative range of values = 67 * EPS, is small (axis 1) > > showProc.time() Time elapsed: 0.26 0.05 0.39 > > 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: 1, Adjusted R-squared: 1 F-statistic: 7.853e+13 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: 1.359e+16 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. 1.332e-15 5.579e-15 9.992e-15 1.623e-14 2.104e-14 6.717e-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.68599e-16" > stopifnot(all.equal(z1[1:10], z1true, tol = 1e-15)) > showProc.time() Time elapsed: 0.19 0.01 0.22 > ##> > ##> 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") Error in axis(side = side, at = at, labels = labels, ...) : CreateAtVector [log-axis()]: axp[0] = 0 < 0! Calls: curve ... plot.default -> localAxis -> Axis -> Axis.default -> axis In addition: Warning messages: 1: In xy.coords(x, y, xlabel, ylabel, log) : 514 y values <= 0 omitted from logarithmic plot 2: In plot.window(...) : Internal(pretty()): very small range.. corrected 3: In axis(side = side, at = at, labels = labels, ...) : CreateAtVector "log"(from axis()): axp[0] = 0 ! Execution halted * checking for unstated dependencies in vignettes ... OK * checking package vignettes in 'inst/doc' ... OK * checking re-building of vignette outputs ... [22s] OK * checking PDF version of manual ... OK * DONE Status: 2 ERRORs, 1 WARNING