R Under development (unstable) (2023-11-02 r85465 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > library(DPQ) > > ### Log-scale {was part of my function code -- should work in R after 2014 :} > lp <- -10*(150:20) > qb <- qbeta(lp, 2,3, log.p=TRUE) > pb <- pbeta(qb, 2,3, log.p=TRUE) > all.equal(lp, pb, countEQ=TRUE) # 0.003525 ... not particularly good; now even *worse*: 0.01925502 [1] "Mean relative difference: 0.01925502" > ## here's the reason : qbeta() stays at minimum > 0, but it should rather underflow > cbind(lp, qb, pb, D = lp-pb, relD = 1 - pb/lp) lp qb pb D relD [1,] -1500 1.112537e-308 -1416.387 -83.612628 0.055741752 [2,] -1490 1.112537e-308 -1416.387 -73.612628 0.049404448 [3,] -1480 1.112537e-308 -1416.387 -63.612628 0.042981505 [4,] -1470 1.112537e-308 -1416.387 -53.612628 0.036471176 [5,] -1460 1.112537e-308 -1416.387 -43.612628 0.029871663 [6,] -1450 1.112537e-308 -1416.387 -33.612628 0.023181123 [7,] -1440 1.112537e-308 -1416.387 -23.612628 0.016397658 [8,] -1430 1.112537e-308 -1416.387 -13.612628 0.009519320 [9,] -1420 1.112537e-308 -1416.387 -3.612628 0.002544104 [10,] -1410 2.712156e-307 -1410.000 0.000000 0.000000000 [11,] -1400 4.025196e-305 -1400.000 0.000000 0.000000000 [12,] -1390 5.973921e-303 -1390.000 0.000000 0.000000000 [13,] -1380 8.866084e-301 -1380.000 0.000000 0.000000000 [14,] -1370 1.315844e-298 -1370.000 0.000000 0.000000000 [15,] -1360 1.952885e-296 -1360.000 0.000000 0.000000000 [16,] -1350 2.898338e-294 -1350.000 0.000000 0.000000000 [17,] -1340 4.301516e-292 -1340.000 0.000000 0.000000000 [18,] -1330 6.384015e-290 -1330.000 0.000000 0.000000000 [19,] -1320 9.474719e-288 -1320.000 0.000000 0.000000000 [20,] -1310 1.406173e-285 -1310.000 0.000000 0.000000000 [21,] -1300 2.086946e-283 -1300.000 0.000000 0.000000000 [22,] -1290 3.097302e-281 -1290.000 0.000000 0.000000000 [23,] -1280 4.596804e-279 -1280.000 0.000000 0.000000000 [24,] -1270 6.822262e-277 -1270.000 0.000000 0.000000000 [25,] -1260 1.012513e-274 -1260.000 0.000000 0.000000000 [26,] -1250 1.502703e-272 -1250.000 0.000000 0.000000000 [27,] -1240 2.230209e-270 -1240.000 0.000000 0.000000000 [28,] -1230 3.309924e-268 -1230.000 0.000000 0.000000000 [29,] -1220 4.912363e-266 -1220.000 0.000000 0.000000000 [30,] -1210 7.290592e-264 -1210.000 0.000000 0.000000000 [31,] -1200 1.082020e-261 -1200.000 0.000000 0.000000000 [32,] -1190 1.605860e-259 -1190.000 0.000000 0.000000000 [33,] -1180 2.383307e-257 -1180.000 0.000000 0.000000000 [34,] -1170 3.537142e-255 -1170.000 0.000000 0.000000000 [35,] -1160 5.249584e-253 -1160.000 0.000000 0.000000000 [36,] -1150 7.791073e-251 -1150.000 0.000000 0.000000000 [37,] -1140 1.156298e-248 -1140.000 0.000000 0.000000000 [38,] -1130 1.716098e-246 -1130.000 0.000000 0.000000000 [39,] -1120 2.546915e-244 -1120.000 0.000000 0.000000000 [40,] -1110 3.779958e-242 -1110.000 0.000000 0.000000000 [41,] -1100 5.609954e-240 -1100.000 0.000000 0.000000000 [42,] -1090 8.325910e-238 -1090.000 0.000000 0.000000000 [43,] -1080 1.235675e-235 -1080.000 0.000000 0.000000000 [44,] -1070 1.833904e-233 -1070.000 0.000000 0.000000000 [45,] -1060 2.721755e-231 -1060.000 0.000000 0.000000000 [46,] -1050 4.039442e-229 -1050.000 0.000000 0.000000000 [47,] -1040 5.995063e-227 -1040.000 0.000000 0.000000000 [48,] -1030 8.897463e-225 -1030.000 0.000000 0.000000000 [49,] -1020 1.320501e-222 -1020.000 0.000000 0.000000000 [50,] -1010 1.959797e-220 -1010.000 0.000000 0.000000000 [51,] -1000 2.908596e-218 -1000.000 0.000000 0.000000000 [52,] -990 4.316739e-216 -990.000 0.000000 0.000000000 [53,] -980 6.406609e-214 -980.000 0.000000 0.000000000 [54,] -970 9.508251e-212 -970.000 0.000000 0.000000000 [55,] -960 1.411150e-209 -960.000 0.000000 0.000000000 [56,] -950 2.094332e-207 -950.000 0.000000 0.000000000 [57,] -940 3.108264e-205 -940.000 0.000000 0.000000000 [58,] -930 4.613073e-203 -930.000 0.000000 0.000000000 [59,] -920 6.846407e-201 -920.000 0.000000 0.000000000 [60,] -910 1.016097e-198 -910.000 0.000000 0.000000000 [61,] -900 1.508021e-196 -900.000 0.000000 0.000000000 [62,] -890 2.238102e-194 -890.000 0.000000 0.000000000 [63,] -880 3.321638e-192 -880.000 0.000000 0.000000000 [64,] -870 4.929748e-190 -870.000 0.000000 0.000000000 [65,] -860 7.316395e-188 -860.000 0.000000 0.000000000 [66,] -850 1.085849e-185 -850.000 0.000000 0.000000000 [67,] -840 1.611543e-183 -840.000 0.000000 0.000000000 [68,] -830 2.391742e-181 -830.000 0.000000 0.000000000 [69,] -820 3.549660e-179 -820.000 0.000000 0.000000000 [70,] -810 5.268163e-177 -810.000 0.000000 0.000000000 [71,] -800 7.818647e-175 -800.000 0.000000 0.000000000 [72,] -790 1.160390e-172 -790.000 0.000000 0.000000000 [73,] -780 1.722172e-170 -780.000 0.000000 0.000000000 [74,] -770 2.555929e-168 -770.000 0.000000 0.000000000 [75,] -760 3.793335e-166 -760.000 0.000000 0.000000000 [76,] -750 5.629809e-164 -750.000 0.000000 0.000000000 [77,] -740 8.355377e-162 -740.000 0.000000 0.000000000 [78,] -730 1.240048e-159 -730.000 0.000000 0.000000000 [79,] -720 1.840394e-157 -720.000 0.000000 0.000000000 [80,] -710 2.731387e-155 -710.000 0.000000 0.000000000 [81,] -700 4.053738e-153 -700.000 0.000000 0.000000000 [82,] -690 6.016281e-151 -690.000 0.000000 0.000000000 [83,] -680 8.928953e-149 -680.000 0.000000 0.000000000 [84,] -670 1.325174e-146 -670.000 0.000000 0.000000000 [85,] -660 1.966733e-144 -660.000 0.000000 0.000000000 [86,] -650 2.918890e-142 -650.000 0.000000 0.000000000 [87,] -640 4.332017e-140 -640.000 0.000000 0.000000000 [88,] -630 6.429283e-138 -630.000 0.000000 0.000000000 [89,] -620 9.541903e-136 -620.000 0.000000 0.000000000 [90,] -610 1.416144e-133 -610.000 0.000000 0.000000000 [91,] -600 2.101744e-131 -600.000 0.000000 0.000000000 [92,] -590 3.119265e-129 -590.000 0.000000 0.000000000 [93,] -580 4.629399e-127 -580.000 0.000000 0.000000000 [94,] -570 6.870637e-125 -570.000 0.000000 0.000000000 [95,] -560 1.019693e-122 -560.000 0.000000 0.000000000 [96,] -550 1.513359e-120 -550.000 0.000000 0.000000000 [97,] -540 2.246023e-118 -540.000 0.000000 0.000000000 [98,] -530 3.333394e-116 -530.000 0.000000 0.000000000 [99,] -520 4.947196e-114 -520.000 0.000000 0.000000000 [100,] -510 7.342289e-112 -510.000 0.000000 0.000000000 [101,] -500 1.089692e-109 -500.000 0.000000 0.000000000 [102,] -490 1.617247e-107 -490.000 0.000000 0.000000000 [103,] -480 2.400207e-105 -480.000 0.000000 0.000000000 [104,] -470 3.562223e-103 -470.000 0.000000 0.000000000 [105,] -460 5.286808e-101 -460.000 0.000000 0.000000000 [106,] -450 7.846319e-99 -450.000 0.000000 0.000000000 [107,] -440 1.164497e-96 -440.000 0.000000 0.000000000 [108,] -430 1.728267e-94 -430.000 0.000000 0.000000000 [109,] -420 2.564975e-92 -420.000 0.000000 0.000000000 [110,] -410 3.806761e-90 -410.000 0.000000 0.000000000 [111,] -400 5.649734e-88 -400.000 0.000000 0.000000000 [112,] -390 8.384949e-86 -390.000 0.000000 0.000000000 [113,] -380 1.244437e-83 -380.000 0.000000 0.000000000 [114,] -370 1.846908e-81 -370.000 0.000000 0.000000000 [115,] -360 2.741054e-79 -360.000 0.000000 0.000000000 [116,] -350 4.068085e-77 -350.000 0.000000 0.000000000 [117,] -340 6.037574e-75 -340.000 0.000000 0.000000000 [118,] -330 8.960554e-73 -330.000 0.000000 0.000000000 [119,] -320 1.329864e-70 -320.000 0.000000 0.000000000 [120,] -310 1.973693e-68 -310.000 0.000000 0.000000000 [121,] -300 2.929221e-66 -300.000 0.000000 0.000000000 [122,] -290 4.347349e-64 -290.000 0.000000 0.000000000 [123,] -280 6.452038e-62 -280.000 0.000000 0.000000000 [124,] -270 9.575673e-60 -270.000 0.000000 0.000000000 [125,] -260 1.421156e-57 -260.000 0.000000 0.000000000 [126,] -250 2.109182e-55 -250.000 0.000000 0.000000000 [127,] -240 3.130304e-53 -240.000 0.000000 0.000000000 [128,] -230 4.645783e-51 -230.000 0.000000 0.000000000 [129,] -220 6.894954e-49 -220.000 0.000000 0.000000000 [130,] -210 1.023302e-46 -210.000 0.000000 0.000000000 [131,] -200 1.518715e-44 -200.000 0.000000 0.000000000 > > require("graphics") > > plot(lp, qb, xlab = quote(log(p)), + main = "qbeta(lp, 2,3, log=TRUE)", log="y", type="l") > ## in the limit, it's fine: > stopifnot(identical(0, qbeta(-Inf, 2,3, log.p=TRUE))) > > ## starting from subnormal > qq <- 2^(-1074:0) > pb. <- pbeta(qq, 2, 3, log.p=TRUE) > ## like a straight line > plot(qq, pb., ylab="pbeta(*, 2, 3, log.p=TRUE)", + log="x", type="l", xaxt="n") > sfsmisc::eaxis(1) > > proc.time() user system elapsed 0.37 0.03 0.39