R Under development (unstable) (2024-08-21 r87038 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > #### History(RCS): /u/maechler/R/MM/NUMERICS/dpq-functions/beta-gamma-etc/qbeta-dist.R,v > > #### Testing qbeta(.), pbeta(.), qt(.), ..... > #### ---------------- > > source(system.file(package="DPQ", "test-tools.R", + mustWork=TRUE))# ../inst/test-tools.R > ## => showProc.time(), .. {from Matrix}, > ## list_() , loadList() , readRDS_() , save2RDS() > > (doExtras <- DPQ:::doExtras()) [1] FALSE > ## save directory (to read from): > (sdir <- system.file("safe", package="DPQ")) [1] "D:/RCompile/CRANincoming/R-devel/lib/DPQ/safe" > > if(!dev.interactive(orNone=TRUE)) pdf("qbeta-dist.pdf") > .O.P. <- par(no.readonly=TRUE) > op <- options(nwarnings = 1e5) > > fcat <- function(..., f.dig= 4, f.wid = f.dig +5, f.flag = ' ', nl = TRUE, + file = "", sep = " ", + fill = FALSE, labels = NULL, append = FALSE) + { + ## Purpose: Formatted CAT -- for printing 'tables' + ## ---------------------------------------------------------------------- + ## Author: Martin Maechler, Date: 12 May 97 + l <- unlist(lapply(list(...), formatC, + wid= f.wid, digits= f.dig, flag= f.flag)) + cat(l, if(nl)"\n", file=file, sep=sep, fill=fill,labels=labels, append=append) + } > > format01prec <- function(x, digits = getOption("digits"), width = digits + 2, + eps = 1e-6, ..., + FUN = function(x,...) formatC(x, flag='-',...)) + { + ## Purpose: format numbers in [0,1] with "precise" result, + ## using "1-.." if necessary. + ## ------------------------------------------------------------------------- + ## Arguments: x: numbers in [0,1]; (still works if not) + ## digits, width: number of digits, width to use with 'FUN' + ## eps: Use '1-' iff x in (1-eps, 1] -- 1e-6 is OPTIMAL + ## ------------------------------------------------------------------------- + ## Author: Martin Maechler, Date: 14 May 97, 18:07 + if(as.integer(digits) < 4) stop('digits must be >= 4') + if(eps < 0 || eps > .1) stop('eps must be in [0, .1]') + i.swap <- is.na(x) | (1-eps < x & x <= 1) #-- Use "1- ." if <~ 1, normal 'FUN' otherwise + r <- character(length(x)) + if(any(i.swap)) + r[i.swap] <- + paste("1-", FUN(1-x[i.swap], digits=digits - 5, width=width-2, ...), + sep='')# -5: '1-' + 4 for exponent -1 for '0' (in other case) + if(any(!i.swap)) + r[!i.swap] <- FUN(x[!i.swap], digits=digits, width=width,...) + attributes(r) <- attributes(x) + r + } > > > ###---------------------- short pre-millennium examples: ---------------- > > ## = ./beta-test-0.R, Jul 9, 1997 (!) > showProc.time() Time (user system elapsed): 0.05 0.03 0.08 > > qbeta(1e-3, 1,10^(1:10)) [1] 1.000450e-04 1.000495e-05 1.000500e-06 1.000500e-07 1.000500e-08 [6] 1.000500e-09 1.000500e-10 1.000500e-11 1.000500e-12 1.000500e-13 > qbeta(1e-5, 1,10^(1:10)) [1] 1.000005e-06 1.000005e-07 1.000005e-08 1.000005e-09 1.000005e-10 [6] 1.000005e-11 1.000005e-12 1.000005e-13 1.000005e-14 1.000005e-15 > qbeta(3e-2, 1,10^(1:10)) [1] 3.041287e-03 3.045457e-04 3.045874e-05 3.045916e-06 3.045920e-07 [6] 3.045921e-08 3.045921e-09 3.045921e-10 3.045921e-11 3.045921e-12 > qbeta(2e-2, 1,10^(1:10)) [1] 2.018231e-03 2.020067e-04 2.020250e-05 2.020269e-06 2.020271e-07 [6] 2.020271e-08 2.020271e-09 2.020271e-10 2.020271e-11 2.020271e-12 > > ##-- quite problematic ones : > cbind(qbeta(2e-2, 1, 10^(0:32))) [,1] [1,] 2.000000e-02 [2,] 2.018231e-03 [3,] 2.020067e-04 [4,] 2.020250e-05 [5,] 2.020269e-06 [6,] 2.020271e-07 [7,] 2.020271e-08 [8,] 2.020271e-09 [9,] 2.020271e-10 [10,] 2.020271e-11 [11,] 2.020271e-12 [12,] 2.020271e-13 [13,] 2.020271e-14 [14,] 2.020271e-15 [15,] 2.020271e-16 [16,] 2.020271e-17 [17,] 2.020271e-18 [18,] 2.020271e-19 [19,] 2.020271e-20 [20,] 2.020271e-21 [21,] 2.020271e-22 [22,] 2.020271e-23 [23,] 2.020271e-24 [24,] 2.020271e-25 [25,] 2.020271e-26 [26,] 2.020271e-27 [27,] 2.020271e-28 [28,] 2.020271e-29 [29,] 2.020271e-30 [30,] 2.020271e-31 [31,] 2.020271e-32 [32,] 2.020271e-33 [33,] 2.020271e-34 > cbind(qbeta(.2, 1, 10^(0:32))) [,1] [1,] 2.000000e-01 [2,] 2.206723e-02 [3,] 2.228948e-03 [4,] 2.231187e-04 [5,] 2.231411e-05 [6,] 2.231433e-06 [7,] 2.231435e-07 [8,] 2.231435e-08 [9,] 2.231436e-09 [10,] 2.231436e-10 [11,] 2.231436e-11 [12,] 2.231436e-12 [13,] 2.231436e-13 [14,] 2.231436e-14 [15,] 2.231436e-15 [16,] 2.231436e-16 [17,] 2.231436e-17 [18,] 2.231436e-18 [19,] 2.231436e-19 [20,] 2.231436e-20 [21,] 2.231436e-21 [22,] 2.231436e-22 [23,] 2.231436e-23 [24,] 2.231436e-24 [25,] 2.231436e-25 [26,] 2.231436e-26 [27,] 2.231436e-27 [28,] 2.231436e-28 [29,] 2.231436e-29 [30,] 2.231436e-30 [31,] 2.231436e-31 [32,] 2.231436e-32 [33,] 2.231436e-33 > cbind(qbeta(.2, .1, 10^(0:32))) [,1] [1,] 1.024000e-07 [2,] 6.508895e-09 [3,] 6.246887e-10 [4,] 6.221601e-11 [5,] 6.219082e-12 [6,] 6.218830e-13 [7,] 6.218805e-14 [8,] 6.218802e-15 [9,] 6.218802e-16 [10,] 6.218802e-17 [11,] 6.218802e-18 [12,] 6.218802e-19 [13,] 6.218802e-20 [14,] 6.218802e-21 [15,] 6.218802e-22 [16,] 6.218802e-23 [17,] 6.218802e-24 [18,] 6.218802e-25 [19,] 6.218802e-26 [20,] 6.218802e-27 [21,] 6.218802e-28 [22,] 6.218802e-29 [23,] 6.218802e-30 [24,] 6.218802e-31 [25,] 6.218802e-32 [26,] 6.218802e-33 [27,] 6.218802e-34 [28,] 6.218802e-35 [29,] 6.218802e-36 [30,] 6.218802e-37 [31,] 6.218802e-38 [32,] 6.218802e-39 [33,] 6.218802e-40 > cbind(qbeta(.1, 1e-1, 10^(0:32)))# "good" [,1] [1,] 1.000000e-10 [2,] 6.356343e-12 [3,] 6.100475e-13 [4,] 6.075782e-14 [5,] 6.073322e-15 [6,] 6.073076e-16 [7,] 6.073051e-17 [8,] 6.073049e-18 [9,] 6.073048e-19 [10,] 6.073048e-20 [11,] 6.073048e-21 [12,] 6.073048e-22 [13,] 6.073048e-23 [14,] 6.073048e-24 [15,] 6.073048e-25 [16,] 6.073048e-26 [17,] 6.073048e-27 [18,] 6.073048e-28 [19,] 6.073048e-29 [20,] 6.073048e-30 [21,] 6.073048e-31 [22,] 6.073048e-32 [23,] 6.073048e-33 [24,] 6.073048e-34 [25,] 6.073048e-35 [26,] 6.073048e-36 [27,] 6.073048e-37 [28,] 6.073048e-38 [29,] 6.073048e-39 [30,] 6.073048e-40 [31,] 6.073048e-41 [32,] 6.073048e-42 [33,] 6.073048e-43 > cbind(qbeta(.1, 1e-2, 10^(0:32)))# look good [,1] [1,] 1.000000e-100 [2,] 5.952797e-102 [3,] 5.688874e-103 [4,] 5.663541e-104 [5,] 5.661018e-105 [6,] 5.660766e-106 [7,] 5.660741e-107 [8,] 5.660738e-108 [9,] 5.660738e-109 [10,] 5.660738e-110 [11,] 5.660738e-111 [12,] 5.660738e-112 [13,] 5.660738e-113 [14,] 5.660738e-114 [15,] 5.660738e-115 [16,] 5.660738e-116 [17,] 5.660738e-117 [18,] 5.660738e-118 [19,] 5.660738e-119 [20,] 5.660738e-120 [21,] 5.660738e-121 [22,] 5.660738e-122 [23,] 5.660738e-123 [24,] 5.660738e-124 [25,] 5.660738e-125 [26,] 5.660738e-126 [27,] 5.660738e-127 [28,] 5.660738e-128 [29,] 5.660738e-129 [30,] 5.660738e-130 [31,] 5.660738e-131 [32,] 5.660738e-132 [33,] 5.660738e-133 > cbind(qbeta(.1, 1e-2, 10^(0:64)))# look good [,1] [1,] 1.000000e-100 [2,] 5.952797e-102 [3,] 5.688874e-103 [4,] 5.663541e-104 [5,] 5.661018e-105 [6,] 5.660766e-106 [7,] 5.660741e-107 [8,] 5.660738e-108 [9,] 5.660738e-109 [10,] 5.660738e-110 [11,] 5.660738e-111 [12,] 5.660738e-112 [13,] 5.660738e-113 [14,] 5.660738e-114 [15,] 5.660738e-115 [16,] 5.660738e-116 [17,] 5.660738e-117 [18,] 5.660738e-118 [19,] 5.660738e-119 [20,] 5.660738e-120 [21,] 5.660738e-121 [22,] 5.660738e-122 [23,] 5.660738e-123 [24,] 5.660738e-124 [25,] 5.660738e-125 [26,] 5.660738e-126 [27,] 5.660738e-127 [28,] 5.660738e-128 [29,] 5.660738e-129 [30,] 5.660738e-130 [31,] 5.660738e-131 [32,] 5.660738e-132 [33,] 5.660738e-133 [34,] 5.660738e-134 [35,] 5.660738e-135 [36,] 5.660738e-136 [37,] 5.660738e-137 [38,] 5.660738e-138 [39,] 5.660738e-139 [40,] 5.660738e-140 [41,] 5.660738e-141 [42,] 5.660738e-142 [43,] 5.660738e-143 [44,] 5.660738e-144 [45,] 5.660738e-145 [46,] 5.660738e-146 [47,] 5.660738e-147 [48,] 5.660738e-148 [49,] 5.660738e-149 [50,] 5.660738e-150 [51,] 5.660738e-151 [52,] 5.660738e-152 [53,] 5.660738e-153 [54,] 5.660738e-154 [55,] 5.660738e-155 [56,] 5.660738e-156 [57,] 5.660738e-157 [58,] 5.660738e-158 [59,] 5.660738e-159 [60,] 5.660738e-160 [61,] 5.660738e-161 [62,] 5.660738e-162 [63,] 5.660738e-163 [64,] 5.660738e-164 [65,] 5.660738e-165 > > plot(qbeta(1/8, 2^-7, 10^(0:64)), log="y")# look good > plot(qbeta(1/8, 2^-8, 10^(0:64)), log="y")# look good (down to 10^-298) > all(qbeta(1/8, 2^-9, 10^(0:64)) == 0)# TRUE, but ... [1] TRUE > all(qbeta(1/8, 2^-10, 10^(0:64)) == 0)# TRUE, but should be positive (in principle) [1] TRUE > > stopifnot(qbeta(.1, 1e-3, 10^(0:32)) == 0) # and we cannot do better: > > curve(pbeta(x, 1e-3, 1 ), 1e-310, log="x", col=2, xaxt="n")# to give nicer axis: > sfsmisc::eaxis(1); axis(1, at=1); abline(h=1,v=1, lty=2) > curve(pbeta(x, 1e-3, 1000), add=TRUE, col=3) > curve(pbeta(x, 1e-3, 1e10), add=TRUE, col=4) > curve(pbeta(x, 1e-3, 1e20), add=TRUE, col=5) > curve(pbeta(x, 1e-3, 1e50), add=TRUE, col=6) > curve(pbeta(x, 1e-3,1e100), add=TRUE, col=5) > curve(pbeta(x, 1e-3,1e200), add=TRUE, col=4)## ==> Warnings already There were 15 warnings (use warnings() to see them) > curve(pbeta(x, 1e-3,1e300), add=TRUE, col=2)## ==> WARNINGS: There were 48 warnings (use warnings() to see them) > ## 1: In pbeta(x, ...) : bgrat(a=1e+300, b=0.001, x=1) *no* convergence: NOTIFY R-core! > summary(warnings()) # ---------------- ^^^^^^^^^^^^^^ Summary of (a total of 48) warning messages: 45x : In pbeta(x, 0.001, 1e+300) : bgrat(a=1e+300, b=0.001, x=1) *no* convergence: NOTIFY R-core! dj=nan, rel.err=nan 1x : In pbeta(x, 0.001, 1e+300) : bgrat(a=1e+300, b=0.001, x=0.999999) *no* convergence: NOTIFY R-core! dj=nan, rel.err=nan 1x : In pbeta(x, 0.001, 1e+300) : bgrat(a=1e+300, b=0.001, x=0.999206) *no* convergence: NOTIFY R-core! dj=nan, rel.err=nan 1x : In pbeta(x, 0.001, 1e+300) : NaNs produced > > ## = ./beta-qbeta-mintst.R, July 31, 1997 (!) > 1 - qbeta(.05, 5e10, 1/2) #-- infinite loop in 0.49; 0.50-a1 [1] 3.84146e-11 > ##-- 3.8416e-11 with MM's qbeta(.) > qbeta(.95, 1/2, 5e10) # the "same" (with swap_tail) [1] 3.841459e-11 > showProc.time() Time (user system elapsed): 0 0.03 0.03 > > ###----------------------- qt(p, df --> Inf) --> qnorm(p) ----------------- > > p.<- 5* 10^-c(6:8,10,15,20,100,300)#- well ... > p <- p. #- try things below with this 'EXTREME' p > p <- c(.975, .995, .9995, 1 - 5e-6) > df1 <- c(1:5,10^c(1:21,25,10*(3:10))) > df0 <- df1[df1< 10^30] ## df0 works for R 0.49 -patch2- > df0f<- unique(sort(c((1:30)/10, df0))) > > > ## system.time() gives 0.0 nowadays ... > sys.time <- function(exp) system.time(for(i in 1:5000) exp) > dig <- 18 > dig <- 12 > dig <- 8 > > cat(" df : CPU[s] | p=", formatC( p , form='e', dig=5, wid=dig+6),"\n") df : CPU[s] | p= 9.75000e-01 9.95000e-01 9.99500e-01 9.99995e-01 > for(df in df0) + cat(formatC(df,w=6),":", formatC(sys.time(qq <- qt(p, df))[1],form='f'), "|", + formatC(qq, form='f', dig=dig,wid=dig+6),"\n") 1 : 0.0000 | 12.70620474 63.65674116 636.61924877 63661.97723111 2 : 0.0000 | 4.30265273 9.92484320 31.59905458 316.22539430 3 : 0.0000 | 3.18244631 5.84090931 12.92397864 60.39682404 4 : 0.0000 | 2.77644511 4.60409487 8.61030158 27.77164766 5 : 0.0000 | 2.57058184 4.03214298 6.86882663 17.89686614 10 : 0.0000 | 2.22813885 3.16927267 4.58689386 8.15028656 100 : 0.0000 | 1.98397152 2.62589052 3.39049131 4.65424006 1000 : 0.0000 | 1.96233908 2.58075470 3.30028265 4.43992647 1e+04 : 0.0000 | 1.96020124 2.57632105 3.29149997 4.41943950 1e+05 : 0.0000 | 1.95998771 2.57587847 3.29062403 4.41739993 1e+06 : 0.0000 | 1.95996636 2.57583422 3.29053646 4.41719606 1e+07 : 0.0000 | 1.95996422 2.57582980 3.29052770 4.41717568 1e+08 : 0.0000 | 1.95996401 2.57582935 3.29052683 4.41717364 1e+09 : 0.0000 | 1.95996399 2.57582931 3.29052674 4.41717344 1e+10 : 0.0000 | 1.95996398 2.57582930 3.29052673 4.41717342 1e+11 : 0.0000 | 1.95996398 2.57582930 3.29052673 4.41717341 1e+12 : 0.0000 | 1.95996398 2.57582930 3.29052673 4.41717341 1e+13 : 0.0000 | 1.95996398 2.57582930 3.29052673 4.41717341 1e+14 : 0.0000 | 1.95996398 2.57582930 3.29052673 4.41717341 1e+15 : 0.0000 | 1.95996398 2.57582930 3.29052673 4.41717341 1e+16 : 0.0000 | 1.95996398 2.57582930 3.29052673 4.41717341 1e+17 : 0.0000 | 1.95996398 2.57582930 3.29052673 4.41717341 1e+18 : 0.0000 | 1.95996398 2.57582930 3.29052673 4.41717341 1e+19 : 0.0000 | 1.95996398 2.57582930 3.29052673 4.41717341 1e+20 : 0.0000 | 1.95996398 2.57582930 3.29052673 4.41717341 1e+21 : 0.0000 | 1.95996398 2.57582930 3.29052673 4.41717341 1e+25 : 0.0000 | 1.95996398 2.57582930 3.29052673 4.41717341 > cat(" _Inf_ : --- |", formatC(qnorm(p), form='f', dig=dig,wid=dig+6),"\n") _Inf_ : --- | 1.95996398 2.57582930 3.29052673 4.41717341 > > > ### Use F-distribution : `` F_{1,nu} == (t_{nu}) ^2 '' > ### Therefore: qt(p,n) == sqrt(qf(2*p-1, 1,n)) for p >= 1/2 > ### ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -- was true for ORIGINAL in R 0.49 > ##for(df in df0f) { ##-- ~ TRUE for (0.49) qt/qf which worked with qbeta > if(any(p < 1/2)) stop("must have p >= 1/2 for these examples!") > > for(df in df0f) { + cat(formatC(df,w = 6),":", all.equal( sqrt(qf(2*p-1, 1,df)), qt(p,df)),"\n") + } 0.1 : TRUE 0.2 : TRUE 0.3 : TRUE 0.4 : TRUE 0.5 : TRUE 0.6 : TRUE 0.7 : TRUE 0.8 : TRUE 0.9 : TRUE 1 : TRUE 1.1 : TRUE 1.2 : TRUE 1.3 : TRUE 1.4 : TRUE 1.5 : TRUE 1.6 : TRUE 1.7 : TRUE 1.8 : TRUE 1.9 : TRUE 2 : TRUE 2.1 : TRUE 2.2 : TRUE 2.3 : TRUE 2.4 : TRUE 2.5 : TRUE 2.6 : TRUE 2.7 : TRUE 2.8 : TRUE 2.9 : TRUE 3 : TRUE 4 : TRUE 5 : TRUE 10 : TRUE 100 : TRUE 1000 : TRUE 1e+04 : TRUE 1e+05 : TRUE 1e+06 : Mean relative difference: 3.240032e-06 1e+07 : Mean relative difference: 3.240022e-07 1e+08 : Mean relative difference: 3.240021e-08 1e+09 : TRUE 1e+10 : TRUE 1e+11 : TRUE 1e+12 : TRUE 1e+13 : TRUE 1e+14 : TRUE 1e+15 : TRUE 1e+16 : TRUE 1e+17 : TRUE 1e+18 : TRUE 1e+19 : TRUE 1e+20 : TRUE 1e+21 : TRUE 1e+25 : TRUE > ## R 3.0.3 : all TRUE, but > ## 1e+06 : Mean relative difference: 3.240032e-06 > ## 1e+07 : Mean relative difference: 3.240022e-07 > ## 1e+08 : Mean relative difference: 3.240021e-08 > ## ==> Applying pt(*) ==> (R 3.5.1) shows that qt() is *better* than qf(.) equivalence > > ## It seems, qt(.) is perfect > p <- 1 - 1/64 # exact rational > df <- 2^seq(10,20, by=1/32) > plot(df, qt(p, df=df)- qnorm(p), log="xy", type="l") > plot(df[-1], -diff(qt(p, df=df)), log="xy", type="l") > plot(df[-(1:2)], diff(qt(p, df=df), diff=2), log="xy", type="l") > all.equal( sqrt(qf(2*p-1, 1,df)), qt(p,df)) ## -> 3.15636e-07 (R 3.0.3, 64bit) [1] "Mean relative difference: 3.156359e-07" > ## looking closely: > plot(df, qf(2*p-1, 1,df) - qt(p,df)^2, type="b") > ## shows clear jump at df = 400'000 > df <- 1000*seq(390, 410, by=1/16) > plot(df, qf(2*p-1, 1,df) - qt(p,df)^2, type="l") > ## ~/R/D/r-devel/R/src/nmath/qf.c : it uses qchisq() for df1 <= df2, df2 > 400000 > ## ==> FIXME: qbeta() is *clearly* better here (than "df = Inf" !) > ## ===== qf() is *clearly* improvable > > showProc.time() Time (user system elapsed): 0.45 0.02 0.47 > > ### Really, the problem is qbeta (.): > ## --------- > ## qt (for x > 1/2) __used__ to be defined as qt(x,n) below, but not longer !! > ## qt(x,n) := sqrt((1 / qbeta( 2*(1-x), n/2, 1/2) - 1) * n) > ##--- with DEBUG, shows already BIG problems (needs *much* CPU time for the large df's > if(FALSE) # is too large and slow + for(p in 2^-(30:2)) { + cat("\n\np = ", formatC(p),"\n===========\n") + cat(sprintf("%6s : %5.3s | %11s %11s | %s\n", + "df","cpu","qbeta","qt", "all.equal(sqrt((1/qb - 1)*.), qt(.))")) + for(df in c(1:5,10^c(1:21,25,10*(3:10)))) { + cpu <- system.time(qb <- qbeta(2 * p, df/2, 1/2))[1] + qt. <- qt(p,df, lower.tail=FALSE) + cat(sprintf("%6g : %5.3g | %11g %11g | %s\n", + df, cpu, qb, qt., all.equal(sqrt(df*(1/qb - 1)), qt.))) + } + } > > > > library(DPQ) Attaching package: 'DPQ' The following object is masked _by_ '.GlobalEnv': format01prec > showProc.time() Time (user system elapsed): 0.02 0 0.02 > > ##----------- Test the 'qnormAppr' function as used in qbeta(.) : > > p <- (0:256)/256 > matplot(p, cbind(qnorm(p), qnormAppr(p)), type="l", ylab="") Warning message: In qnormAppr(p) : 'qnormAppr' is deprecated. Use 'qnormUappr' instead. See help("Deprecated") > legend("topleft", c("qnorm(p)", "qnormAppr(p)"), col=1:2, lty=1:2, bty="n") > abline(h=0,v=(0:2)/2, lty=2, col="gray60") > ## absolute error: > curve(qnorm(x) - qnormAppr(x), col=2, lwd=2, n=1001, 0, 1) Warning message: In qnormAppr(x) : 'qnormAppr' is deprecated. Use 'qnormUappr' instead. See help("Deprecated") > abline(h=0,v=(0:2)/2, lty=2, col="gray60") > > ##==> Really make use of symmetry --> use qnormUappr(): > matplot(p, cbind(qnorm(p), qnormUappr(p)), type="l", ylab="") > abline(h=0,v=(0:2)/2, lty=2, col="gray60") > > ## absolute error: > c.1 <- rgb(1,0,0); c.2 <- adjustcolor("black",.5) > curve(qnorm(x) - qnormUappr(x), col=c.1, lwd=2, n=1001, 0, 1) > curve(qnorm(x) - qnormAppr (x), col=c.2, lwd=4, n=1001, + add=TRUE, lty=3) Warning message: In qnormAppr(x) : 'qnormAppr' is deprecated. Use 'qnormUappr' instead. See help("Deprecated") > legend(.1, .0025, paste("qnorm(.) -",c("qnormUappr(.)", " qnormAppr (.)")), + col=c(c.1,c.2), lwd=c(2,3), lty=c(1,3), bty="n") > abline(h=0,v=(0:2)/2, lty=2, col="gray60") > showProc.time() Time (user system elapsed): 0.02 0 0.01 > > > ## From R's source /tests/d-p-q-r-tests.R : > options(rErr.eps = 1e-30) # {*not* to be reset via options(op)} > rErr <- function(approx, true, eps = .Options$rErr.eps) + { + if(is.null(eps)) { eps <- 1e-30; options(rErr.eps = eps) } + ifelse(Mod(true) >= eps, + 1 - approx / true, # relative error + true - approx) # absolute error (e.g. when true=0) + } > > ## The relative error of this approximation is quite Asymmetric: mainly < 0 > x <- c(10^(-9:-4),seq(1e-3, .6, len = 1+2^11)) > plot(1-x, rErr(qnormAppr(1-x), qnorm(1-x)), type = 'l', col = 'red', + main = "Rel. Error of 'qnormAppr(1-x)'") Warning message: In qnormAppr(1 - x) : 'qnormAppr' is deprecated. Use 'qnormUappr' instead. See help("Deprecated") > abline(h = 0, col = 'gray') > ## Note: qnorm(.) used to use AS 111 (26)1977, but has been changed to > ## ---- AS 241 (37)1988 which is more accurate {and a bit slower} > ## > ## Much more sensical in modern R {but looks identical, for 1-x > 1/2 }: > lines(1-x, rErr(qnormUappr(x, lower=FALSE), qnorm(x, lower=FALSE)), + lwd=3, col = adjustcolor("skyblue", 1/2)) > > curve(qnorm(x) - qnormUappr(x), col=2, lwd=2, n=1001, .44, 1) > abline(h=0, col=adjustcolor(1,.5)); axis(1, at=.44) > ## Warning message: > ## In qnormUappr(x) : p[.] < 1/2 & lower tail is inaccurate > curve(qnorm(x) - qnormAppr(x), lwd=5, n=1001, col=adjustcolor(3,1/3), add=TRUE) Warning message: In qnormAppr(x) : 'qnormAppr' is deprecated. Use 'qnormUappr' instead. See help("Deprecated") > ##par(new=TRUE) > ## 1/10 * rel.error(): > cblu <- adjustcolor("blue4",.5) > curve(rErr(qnormUappr(x), qnorm(x)) / 10, n=1001, + col = cblu, lwd=3, lty=2, add=TRUE) > ## axis(4, col=cblu, col.axis=cblu) > showProc.time() Time (user system elapsed): 0.03 0 0.03 > > > ###-------- Testing my own qbeta.R(.) [[and qbetaAppr(.) used there] > > qbeta (.05 , 10/2, 1/2) #-== 0.6682436 [1] 0.6682436 > qbeta (.05 , 100/2, 1/2) #-== 0.9621292 [1] 0.9621292 > qbeta (.05 ,1000/2, 1/2) #-== 0.996164 [1] 0.996164 > > ###=========== FIXME: More qbetaAppr*() versions > ###============================================= > ## and use lower.tail=TRUE/FALSE *and* log.p=FALSE/TRUE > ## and determine regions where qbetaAppr*() are very accurate ((so, say 1 Newton step then is perfect!)) > > qb.names <- paste0("qbeta", c('', '.R', 'Appr', 'Appr.1')) > qf <- lapply(setNames(,qb.names), get) > str(qf) ## for (qn in qb.names) cat(qn,":", deparse(args(get(qn))),"\n\n") List of 4 $ qbeta :function (p, shape1, shape2, ncp = 0, lower.tail = TRUE, log.p = FALSE) $ qbeta.R :function (alpha, p, q, lower.tail = TRUE, log.p = FALSE, logbeta = lbeta(p, q), low.bnd = 3e-308, up.bnd = 1 - 2.22e-16, method = c("AS109", "Newton-log"), tol.outer = 1e-15, f.acu = function(a, p, q) max(1e-300, 10^(-13 - 2.5/pp^2 - 0.5/a^2)), fpu = .Machine$double.xmin, qnormU.fun = function(u, lu) qnormUappr(p = u, lp = lu), R.pre.2014 = FALSE, verbose = getOption("verbose"), non.finite.report = verbose) $ qbetaAppr :function (a, p, q, lower.tail = TRUE, log.p = FALSE, y = qnormUappr(a, lower.tail = lower.tail, log.p = log.p), logbeta = lbeta(p, q), verbose = getOption("verbose") && length(a) == 1) $ qbetaAppr.1:function (a, p, q, lower.tail = TRUE, log.p = FALSE, y = qnormUappr(a, lower.tail = lower.tail, log.p = log.p)) > > p.set2 <- c(1:5,10,20,50,100,500,1000, 10^c(4:20,25,10*3:10)) > p.set2 <- c(1:5,10,20,50,100,500,1000, 10^c(4:18,20,30,100)) > ## TODO: still too expensive > p.set2 <- c(1e-10,1e-5,1e-3,.1,.3,.5, .8, 1,10, 10^c(6:18,20,30,100)) > q.set2 <- c(1e-4, 1/2, 1, 3/2, 2, 5, 1e4) > > pn2 <- paste("p=",formatC(p.set2,wid = 1),sep = '') > qn2 <- paste("q=",formatC(q.set2,wid = 1),sep = '') > > sfil1 <- file.path(sdir, "tests_qbeta-d-ssR.rds") > if(!doExtras && file.exists(sfil1)) { + ssR_l <- readRDS_(sfil1) + str(ssR_l) + loadList(ssR_l) + + } else { ## do run the simulation [always if(doExtras)] : + + ra <- array(dim = c(length(q.set2), length(qb.names), length(p.set2)), + dimnames = list(qn2, qb.names, pn2)) + ind_ct <- names(system.time(1))[ c(1,3) ] # "user" and "elapsed" + cta <- array(dim = c(length(ind_ct), dim(ra)), + dimnames = c(list(ind_ct), dimnames(ra))) + for(qq in q.set2) { + cat("\nq=",qq,"\n======\n") + for (p in p.set2) { + cat(sprintf("p=%-7g : >> ",p)) + for (qn in qb.names) { + cat(" ",qn,": ", sep = '') + st <- system.time( + r <- qf[[qn]](.05, p, qq))[ind_ct] + cta[, qn2[qq == q.set2], qn, pn2[p == p.set2]] <- st + ra [ qn2[qq == q.set2], qn, pn2[p == p.set2]] <- r + }; cat("\n") + } + showProc.time() ## -- similar for all q + } + + print(summary(warnings())) # many warnings from qbeta() inaccuracies + + save2RDS(list_(ra, cta), file = sfil1) + + }## {run simulations} ------------------------------------------------------- Reading from D:/RCompile/CRANincoming/R-devel/lib/DPQ/safe/tests_qbeta-d-ssR.rds Time (user system elapsed): 0 0 0 List of 2 $ ra : num [1:7, 1:4, 1:25] 0 0 0 0 0 ... ..- attr(*, "dimnames")=List of 3 .. ..$ : chr [1:7] "q=0.0001" "q=0.5" "q=1" "q=1.5" ... .. ..$ : chr [1:4] "qbeta" "qbeta.R" "qbetaAppr" "qbetaAppr.1" .. ..$ : chr [1:25] "p=1e-10" "p=1e-05" "p=0.001" "p=0.1" ... $ cta: num [1:2, 1:7, 1:4, 1:25] 0 0 0 0 0 0 0 0 0 0 ... ..- attr(*, "dimnames")=List of 4 .. ..$ : chr [1:2] "user.self" "elapsed" .. ..$ : chr [1:7] "q=0.0001" "q=0.5" "q=1" "q=1.5" ... .. ..$ : chr [1:4] "qbeta" "qbeta.R" "qbetaAppr" "qbetaAppr.1" .. ..$ : chr [1:25] "p=1e-10" "p=1e-05" "p=0.001" "p=0.1" ... > > ## Timings : > 1e4*apply(cta, 2:3, mean) qbeta qbeta.R qbetaAppr qbetaAppr.1 q=0.0001 35.2 18.0 0.6 1.6 q=0.5 36.0 6.4 0.4 0.6 q=1 35.2 6.4 0.8 0.4 q=1.5 35.0 7.4 0.2 0.4 q=2 35.2 6.4 0.4 0.4 q=5 35.2 7.2 0.6 0.2 q=1e+04 35.2 6.8 0.2 0.2 > 1e4*apply(cta, c(2,4), mean) p=1e-10 p=1e-05 p=0.001 p=0.1 p=0.3 p=0.5 p=0.8 p=1 p=10 p=1e+06 q=0.0001 7.50 15.00 82.50 1.25 2.50 1.25 1.25 1.25 1.25 1.25 q=0.5 7.50 12.50 15.00 3.75 0.00 0.00 1.25 0.00 0.00 0.00 q=1 6.25 13.75 15.00 2.50 1.25 0.00 0.00 0.00 0.00 0.00 q=1.5 6.25 13.75 16.25 0.00 2.50 1.25 0.00 1.25 0.00 0.00 q=2 6.25 13.75 16.25 0.00 1.25 0.00 0.00 1.25 0.00 0.00 q=5 7.50 12.50 16.25 0.00 1.25 0.00 0.00 0.00 2.50 0.00 q=1e+04 7.50 13.75 15.00 0.00 1.25 1.25 2.50 0.00 0.00 0.00 p=1e+07 p=1e+08 p=1e+09 p=1e+10 p=1e+11 p=1e+12 p=1e+13 p=1e+14 q=0.0001 1.25 1.25 0.00 0.00 2.50 0.00 2.50 0 q=0.5 0.00 1.25 1.25 1.25 0.00 1.25 1.25 0 q=1 0.00 3.75 1.25 1.25 0.00 0.00 0.00 0 q=1.5 1.25 1.25 1.25 0.00 1.25 1.25 0.00 0 q=2 0.00 0.00 2.50 1.25 0.00 0.00 1.25 0 q=5 0.00 0.00 1.25 0.00 1.25 1.25 0.00 0 q=1e+04 0.00 0.00 0.00 1.25 0.00 1.25 0.00 0 p=1e+15 p=1e+16 p=1e+17 p=1e+18 p=1e+20 p=1e+30 p=1e+100 q=0.0001 1.25 0.00 1.25 2.50 72.50 72.50 73.75 q=0.5 0.00 1.25 0.00 1.25 77.50 72.50 72.50 q=1 0.00 0.00 2.50 0.00 73.75 73.75 72.50 q=1.5 0.00 0.00 0.00 0.00 73.75 73.75 73.75 q=2 0.00 1.25 0.00 1.25 73.75 72.50 72.50 q=5 0.00 5.00 1.25 0.00 75.00 72.50 72.50 q=1e+04 0.00 0.00 1.25 0.00 73.75 72.50 73.75 > > noquote(format01prec(aperm(ra))) , , q=0.0001 qbeta qbeta.R qbetaAppr qbetaAppr.1 p=1e-10 0 3.329564e-303 -Inf 1-NaN p=1e-05 0 3.426311e-301 -Inf 1-NaN p=0.001 2.30472e-260 2.30472e-260 -1.452749e+191 1-NaN p=0.1 1-1.1e-16 1-6.7e-16 1-0 1-NaN p=0.3 1-1.1e-16 1-6.7e-16 1-0 1-NaN p=0.5 1-1.1e-16 1-1e-15 1-0 1-NaN p=0.8 1-1.1e-16 1-1.6e-15 1-0 0.9838125 p=1 1-1.1e-16 1-1.1e-15 1-0 1-NaN p=10 1-1.1e-16 1-7.8e-16 1-0 1-NaN p=1e+06 1-1.1e-16 0.5 1-0 1-NaN p=1e+07 1-1.1e-16 0.5 1-0 1-NaN p=1e+08 1-1.1e-16 0.5 1-0 1-NaN p=1e+09 1-1.1e-16 0.5 1-0 1-NaN p=1e+10 1-1.1e-16 0.5 1-0 1-NaN p=1e+11 1-1.1e-16 0.5 1-0 1-NaN p=1e+12 1-1.1e-16 0.5 1-0 1-NaN p=1e+13 1-1.1e-16 0.5 1-0 1-NaN p=1e+14 1-1.1e-16 0.5 1-0 1-NaN p=1e+15 1-1.1e-16 0.5 1-0 1-NaN p=1e+16 1-1.1e-16 0.5 1-0 1-NaN p=1e+17 1-1.1e-16 0.5 1-0 1-NaN p=1e+18 1-1.1e-16 0.5 1-0 1-NaN p=1e+20 5.562685e-309 0.5 1-0 1-NaN p=1e+30 5.562685e-309 0.5 1-0 1-NaN p=1e+100 5.562685e-309 0.5 1-0 1-NaN , , q=0.5 qbeta qbeta.R qbetaAppr qbetaAppr.1 p=1e-10 0 6.439797e-303 1-0 1-NaN p=1e-05 0 2.345488e-301 1-0 1-NaN p=0.001 0 5.574727e-301 1-0 1-NaN p=0.1 3.383362e-13 3.383362e-13 0.5661543 1-NaN p=0.3 0.0001303367 0.0001303367 1-1.7e-07 1-NaN p=0.5 0.00615583 0.00615583 1-3.3e-08 1-NaN p=0.8 0.04993758 0.04993758 1-1.5e-08 1-NaN p=1 0.0975 0.0975 1-1.1e-08 1-NaN p=10 0.8213133 0.8213133 1-8.6e-10 1-NaN p=1e+06 0.9999981 0.9999981 1-8.3e-15 1-NaN p=1e+07 1-1.9e-07 1-1.9e-07 1-8.9e-16 1-NaN p=1e+08 1-1.9e-08 0.5 1-1.1e-16 1-NaN p=1e+09 1-1.9e-09 0.5 1-0 1-NaN p=1e+10 1-1.9e-10 0.5 1-0 1-NaN p=1e+11 1-1.9e-11 0.5 1-0 1-NaN p=1e+12 1-1.9e-12 0.5 1-0 1-NaN p=1e+13 1-1.9e-13 0.5 1-0 1-NaN p=1e+14 1-1.9e-14 0.5 1-0 1-NaN p=1e+15 1-1.9e-15 0.5 1-0 1-NaN p=1e+16 1-2.2e-16 0.5 1-0 1-NaN p=1e+17 1-1.1e-16 0.5 1-0 1-NaN p=1e+18 1-1.1e-16 0.5 1-0 1-NaN p=1e+20 5.562685e-309 0.5 1-0 1-NaN p=1e+30 5.562685e-309 0.5 1-0 1-NaN p=1e+100 5.562685e-309 0.5 1-0 1-NaN , , q=1 qbeta qbeta.R qbetaAppr qbetaAppr.1 p=1e-10 0 7.914119e-303 1-0 1-NaN p=1e-05 0 2.591748e-301 1-0 1-NaN p=0.001 0 6.29134e-301 1-0 1-NaN p=0.1 9.765625e-14 9.765625e-14 0.6697565 1-NaN p=0.3 4.605039e-05 4.605039e-05 0.876302 1-NaN p=0.5 0.0025 0.0025 0.9238985 1-NaN p=0.8 0.02364354 0.02364354 0.951748 0.0228813 p=1 0.05 0.05 0.9612113 0.03890576 p=10 0.7411344 0.7411344 0.9960522 0.6831777 p=1e+06 0.999997 0.999997 1-4e-08 0.999996 p=1e+07 1-3e-07 1-3e-07 1-4e-09 1-4e-07 p=1e+08 1-3e-08 1-3e-08 1-4e-10 1-4e-08 p=1e+09 1-3e-09 1-3e-09 1-4e-11 1-4e-09 p=1e+10 1-3e-10 1-3e-10 1-4e-12 1-4e-10 p=1e+11 1-3e-11 1-3e-11 1-4e-13 1-4e-11 p=1e+12 1-3e-12 1-3e-12 1-4e-14 1-4e-12 p=1e+13 1-3e-13 1-3e-13 1-4e-15 1-4e-13 p=1e+14 1-3e-14 1-2.9e-14 1-4.4e-16 1-4e-14 p=1e+15 1-3e-15 0.5 1-0 1-4e-15 p=1e+16 1-3.3e-16 0.5 1-0 1-4.4e-16 p=1e+17 1-1.1e-16 0.5 1-0 1-0 p=1e+18 1-1.1e-16 0.5 1-0 1-0 p=1e+20 5.562685e-309 0.5 1-0 1-0 p=1e+30 5.562685e-309 0.5 1-0 1-0 p=1e+100 5.562685e-309 0.5 1-0 1-0 , , q=1.5 qbeta qbeta.R qbetaAppr qbetaAppr.1 p=1e-10 0 6.855128e-303 0.5056046 1-NaN p=1e-05 0 2.523787e-301 0.5056195 1-NaN p=0.001 0 6.07494e-301 0.5070889 1-NaN p=0.1 5.464319e-14 5.464319e-14 0.6200235 1-NaN p=0.3 2.720778e-05 2.720778e-05 0.7402516 1-NaN p=0.5 0.001542919 0.001542919 0.8026842 1-NaN p=0.8 0.01540694 0.01540694 0.8549721 0.01624554 p=1 0.03361747 0.03361747 0.8767465 0.0297495 p=10 0.6827063 0.6827063 0.6718925 0.6718925 p=1e+06 0.9999961 0.9999961 0.9999959 0.9999959 p=1e+07 1-3.9e-07 1-3.9e-07 1-4.1e-07 1-4.1e-07 p=1e+08 1-3.9e-08 1-3.9e-08 1-4.1e-08 1-4.1e-08 p=1e+09 1-3.9e-09 1-3.9e-09 1-4.1e-09 1-4.1e-09 p=1e+10 1-3.9e-10 1-3.9e-10 1-4.1e-10 1-4.1e-10 p=1e+11 1-3.9e-11 1-3.9e-11 1-4.1e-11 1-4.1e-11 p=1e+12 1-3.9e-12 1-3.9e-12 1-4.1e-12 1-4.1e-12 p=1e+13 1-3.9e-13 1-3.9e-13 1-4.1e-13 1-4.1e-13 p=1e+14 1-3.9e-14 1-3.9e-14 1-4.1e-14 1-4.1e-14 p=1e+15 1-3.9e-15 1-4.1e-15 1-4.1e-15 1-4.1e-15 p=1e+16 1-3.3e-16 1-4.4e-16 1-4.4e-16 1-4.4e-16 p=1e+17 1-1.1e-16 0.5 1-0 1-0 p=1e+18 1-1.1e-16 0.5 1-0 1-0 p=1e+20 5.562685e-309 0.5 1-0 1-0 p=1e+30 5.562685e-309 0.5 1-0 1-0 p=1e+100 5.562685e-309 0.5 1-0 1-0 , , q=2 qbeta qbeta.R qbetaAppr qbetaAppr.1 p=1e-10 0 1.738045e-303 0.4866481 1-NaN p=1e-05 0 3.039556e-301 0.4866558 1-NaN p=0.001 0 7.354448e-301 0.4874102 1-NaN p=0.1 3.765071e-14 3.765071e-14 0.5530882 1-NaN p=0.3 1.920563e-05 1.920563e-05 0.6449835 1-NaN p=0.5 0.001111935 0.001111935 0.7055327 1-NaN p=0.8 0.01141254 0.01141254 0.7655198 0.0117015 p=1 0.02532057 0.02532057 0.7935567 0.02216264 p=10 0.6356405 0.6356405 0.631679 0.631679 p=1e+06 0.9999953 0.9999953 0.9999951 0.9999951 p=1e+07 1-4.7e-07 1-4.7e-07 1-4.9e-07 1-4.9e-07 p=1e+08 1-4.7e-08 1-4.7e-08 1-4.9e-08 1-4.9e-08 p=1e+09 1-4.7e-09 1-4.7e-09 1-4.9e-09 1-4.9e-09 p=1e+10 1-4.7e-10 1-4.7e-10 1-4.9e-10 1-4.9e-10 p=1e+11 1-4.7e-11 1-4.7e-11 1-4.9e-11 1-4.9e-11 p=1e+12 1-4.7e-12 1-4.7e-12 1-4.9e-12 1-4.9e-12 p=1e+13 1-4.7e-13 1-4.7e-13 1-4.9e-13 1-4.9e-13 p=1e+14 1-4.7e-14 1-4.7e-14 1-4.9e-14 1-4.9e-14 p=1e+15 1-4.8e-15 1-4.9e-15 1-4.9e-15 1-4.9e-15 p=1e+16 1-4.4e-16 1-4.4e-16 1-4.4e-16 1-4.4e-16 p=1e+17 1-1.1e-16 0.5 1-0 1-0 p=1e+18 1-1.1e-16 0.5 1-0 1-0 p=1e+20 5.562685e-309 0.5 1-0 1-0 p=1e+30 5.562685e-309 0.5 1-0 1-0 p=1e+100 5.562685e-309 0.5 1-0 1-0 , , q=5 qbeta qbeta.R qbetaAppr qbetaAppr.1 p=1e-10 0 8.947603e-304 0.3408882 1-NaN p=1e-05 0 3.714058e-301 0.3408904 1-NaN p=0.001 0 8.115442e-301 0.3411091 1-NaN p=0.1 1.30088e-14 1.30088e-14 0.3622665 1-NaN p=0.3 6.89342e-06 6.89342e-06 0.401116 1-NaN p=0.5 0.0004132543 0.0004132543 0.435504 1-NaN p=0.8 0.004456116 0.004456116 0.4802686 0.004233553 p=1 0.01020622 0.01020622 0.5063654 0.008471312 p=10 0.4599946 0.4599946 0.4599072 0.4599072 p=1e+06 0.9999908 0.9999908 0.9999908 0.9999908 p=1e+07 1-9.2e-07 1-9.2e-07 1-9.2e-07 1-9.2e-07 p=1e+08 1-9.2e-08 1-9.2e-08 1-9.2e-08 1-9.2e-08 p=1e+09 1-9.2e-09 1-9.2e-09 1-9.2e-09 1-9.2e-09 p=1e+10 1-9.2e-10 1-9.2e-10 1-9.2e-10 1-9.2e-10 p=1e+11 1-9.2e-11 1-9.2e-11 1-9.2e-11 1-9.2e-11 p=1e+12 1-9.2e-12 1-9.2e-12 1-9.2e-12 1-9.2e-12 p=1e+13 1-9.2e-13 1-9.2e-13 1-9.2e-13 1-9.2e-13 p=1e+14 1-9.1e-14 1-9.2e-14 1-9.2e-14 1-9.2e-14 p=1e+15 1-9.1e-15 1-9.1e-15 1-9.1e-15 1-9.1e-15 p=1e+16 1-8.9e-16 1-1e-15 1-1e-15 1-1e-15 p=1e+17 1-1.1e-16 0.5 1-1.1e-16 1-1.1e-16 p=1e+18 1-1.1e-16 0.5 1-0 1-0 p=1e+20 5.562685e-309 0.5 1-0 1-0 p=1e+30 5.562685e-309 0.5 1-0 1-0 p=1e+100 5.562685e-309 0.5 1-0 1-0 , , q=1e+04 qbeta qbeta.R qbetaAppr qbetaAppr.1 p=1e-10 0 3.040562e-303 0.008211652 1-NaN p=1e-05 0 2.651973e-301 0.008211653 1-NaN p=0.001 0 5.644618e-301 0.008211752 1-NaN p=0.1 5.930978e-18 5.930978e-18 0.008221652 1-NaN p=0.3 3.211147e-09 3.211147e-09 0.008241652 1-NaN p=0.5 1.966119e-07 1.966119e-07 0.008261651 1-NaN p=0.8 2.189773e-06 2.189773e-06 0.008291648 1.974375e-06 p=1 5.129316e-06 5.129316e-06 0.008311645 4.083281e-06 p=10 0.0005421496 0.0005421496 0.0005408577 0.0005408577 p=1e+06 0.9899364 0.9899364 0.9899364 0.9899364 p=1e+07 0.9989845 0.9989845 0.9989845 0.9989845 p=1e+08 0.9998984 0.9998984 0.9998984 0.9998984 p=1e+09 0.9999898 0.9999898 0.9999898 0.9999898 p=1e+10 0.999999 0.999999 0.999999 0.999999 p=1e+11 1-1e-07 1-1e-07 1-1e-07 1-1e-07 p=1e+12 1-1e-08 1-1e-08 1-1e-08 1-1e-08 p=1e+13 1-1e-09 1-1e-09 1-1e-09 1-1e-09 p=1e+14 1-1e-10 1-1e-10 1-1e-10 1-1e-10 p=1e+15 1-1e-11 1-1e-11 1-1e-11 1-1e-11 p=1e+16 1-1e-12 1-1e-12 1-1e-12 1-1e-12 p=1e+17 1-1e-13 1-1e-13 1-1e-13 1-1e-13 p=1e+18 1-1e-14 1-1e-14 1-1e-14 1-1e-14 p=1e+20 5.562685e-309 0.5 1-1.1e-16 1-1.1e-16 p=1e+30 5.562685e-309 0.5 1-0 1-0 p=1e+100 5.562685e-309 0.5 1-0 1-0 > > qbeta.R(.05 , 10/2, 1/2, verbose=2) no swap logbeta= -0.2073951943460705938 ; acu= 7.94328e-214 p or q <= 1 : t > 0; t' > 1 : qbetaAppr.4() = 1-2/(t'+1) initial 'xinbta' = 0.999999998240895471 y=pbeta(x..)= 0.999896784146897 y.n=(y-a)*e^..= 3.238e-05 Inner loop: < 1 it --> adj=g*y= 3.238e-05 --Outer Loop: |tx - xinbta| > eps; tx= 0.99996762011888396 tx-xinbta = -3.238e-05 y=pbeta(x..)= 0.985997052538556 y.n=(y-a)*e^..= 0.004329 Inner loop: < 1 it --> adj=g*y= 0.004329 --Outer Loop: |tx - xinbta| > eps; tx= 0.995638518621397739 tx-xinbta = -0.004329 y=pbeta(x..)= 0.838417390393245 y.n=(y-a)*e^..= 0.04306 Inner loop: < 1 it --> adj=g*y= 0.04306 --Outer Loop: |tx - xinbta| > eps; tx= 0.952576353883728744 tx-xinbta = -0.04306 y=pbeta(x..)= 0.496555273537214 y.n=(y-a)*e^..= 0.09598 Inner loop: < 1 it --> adj=g*y= 0.09598 --Outer Loop: |tx - xinbta| > eps; tx= 0.856591453633085864 tx-xinbta = -0.09598 y=pbeta(x..)= 0.224784573477933 y.n=(y-a)*e^..= 0.09991 Inner loop: < 1 it --> adj=g*y= 0.09991 --Outer Loop: |tx - xinbta| > eps; tx= 0.756677745075267971 tx-xinbta = -0.09991 y=pbeta(x..)= 0.103183739893054 y.n=(y-a)*e^..= 0.06504 Inner loop: < 1 it --> adj=g*y= 0.06504 --Outer Loop: |tx - xinbta| > eps; tx= 0.691641515161518616 tx-xinbta = -0.06504 y=pbeta(x..)= 0.0608901635936472 y.n=(y-a)*e^..= 0.02148 Inner loop: < 1 it --> adj=g*y= 0.02148 --Outer Loop: |tx - xinbta| > eps; tx= 0.670164815056873309 tx-xinbta = -0.02148 y=pbeta(x..)= 0.0508243370848005 y.n=(y-a)*e^..= 0.001907 Inner loop: < 1 it --> adj=g*y= 0.001907 --Outer Loop: |tx - xinbta| > eps; tx= 0.668257351499830343 tx-xinbta = -0.001907 y=pbeta(x..)= 0.0500058585993167 y.n=(y-a)*e^..= 1.375e-05 Inner loop: < 1 it --> adj=g*y= 1.375e-05 --Outer Loop: |tx - xinbta| > eps; tx= 0.668243600037566221 tx-xinbta = -1.375e-05 y=pbeta(x..)= 0.050000000301821 y.n=(y-a)*e^..= 7.085e-10 Inner loop: < 1 it --> adj=g*y= 7.085e-10 --Outer Loop: |tx - xinbta| > eps; tx= 0.668243599329050864 tx-xinbta = -7.085e-10 y=pbeta(x..)= 0.05 y.n=(y-a)*e^..= 8.144e-17 Inner loop: < 1 it --> adj=g*y= 8.144e-17 10 outer iterations [1] 0.6682436 > qbeta.R(.05 , 100/2, 1/2, verbose=2) no swap logbeta= -1.381146601451041178 ; acu= 9.977e-214 p or q <= 1 : t > 0; t' > 1 : qbetaAppr.4() = 1-2/(t'+1) initial 'xinbta' = 0.999999999832045239 y=pbeta(x..)= 0.999896854448551 y.n=(y-a)*e^..= 3.093e-06 Inner loop: < 1 it --> adj=g*y= 3.093e-06 --Outer Loop: |tx - xinbta| > eps; tx= 0.999996906345258951 tx-xinbta = -3.093e-06 y=pbeta(x..)= 0.986001924843581 y.n=(y-a)*e^..= 0.0004138 Inner loop: < 1 it --> adj=g*y= 0.0004138 --Outer Loop: |tx - xinbta| > eps; tx= 0.999583140964995565 tx-xinbta = -0.0004138 y=pbeta(x..)= 0.838601442528779 y.n=(y-a)*e^..= 0.00413 Inner loop: < 1 it --> adj=g*y= 0.00413 --Outer Loop: |tx - xinbta| > eps; tx= 0.995453611039165853 tx-xinbta = -0.00413 y=pbeta(x..)= 0.50072299189132 y.n=(y-a)*e^..= 0.009547 Inner loop: < 1 it --> adj=g*y= 0.009547 --Outer Loop: |tx - xinbta| > eps; tx= 0.985906139603156206 tx-xinbta = -0.009547 y=pbeta(x..)= 0.234668101887079 y.n=(y-a)*e^..= 0.01104 Inner loop: < 1 it --> adj=g*y= 0.01104 --Outer Loop: |tx - xinbta| > eps; tx= 0.974861805535010451 tx-xinbta = -0.01104 y=pbeta(x..)= 0.111467953610859 y.n=(y-a)*e^..= 0.008527 Inner loop: < 1 it --> adj=g*y= 0.008527 --Outer Loop: |tx - xinbta| > eps; tx= 0.966335090752517889 tx-xinbta = -0.008527 y=pbeta(x..)= 0.0649033923227893 y.n=(y-a)*e^..= 0.00368 Inner loop: < 1 it --> adj=g*y= 0.00368 --Outer Loop: |tx - xinbta| > eps; tx= 0.962655582868875981 tx-xinbta = -0.00368 y=pbeta(x..)= 0.0516510205056109 y.n=(y-a)*e^..= 0.0005176 Inner loop: < 1 it --> adj=g*y= 0.0005176 --Outer Loop: |tx - xinbta| > eps; tx= 0.96213801750618333 tx-xinbta = -0.0005176 y=pbeta(x..)= 0.0500271442087208 y.n=(y-a)*e^..= 8.797e-06 Inner loop: < 1 it --> adj=g*y= 8.797e-06 --Outer Loop: |tx - xinbta| > eps; tx= 0.962129220738687008 tx-xinbta = -8.797e-06 y=pbeta(x..)= 0.0500000076554689 y.n=(y-a)*e^..= 2.482e-09 Inner loop: < 1 it --> adj=g*y= 2.482e-09 --Outer Loop: |tx - xinbta| > eps; tx= 0.96212921825633857 tx-xinbta = -2.482e-09 y=pbeta(x..)= 0.0500000000000005 y.n=(y-a)*e^..= 1.597e-16 Inner loop: < 1 it --> adj=g*y= 1.597e-16 10 outer iterations [1] 0.9621292 > qbeta.R(.05 ,1000/2, 1/2, verbose=2) no swap logbeta= -2.534689106328062547 ; acu= 9.99977e-214 p or q <= 1 : t > 0; t' > 1 : qbetaAppr.4() = 1-2/(t'+1) initial 'xinbta' = 0.999999999983280152 y=pbeta(x..)= 0.999896855126961 y.n=(y-a)*e^..= 3.08e-07 Inner loop: < 1 it --> adj=g*y= 3.08e-07 --Outer Loop: |tx - xinbta| > eps; tx= 0.99999969202554273 tx-xinbta = -3.08e-07 y=pbeta(x..)= 0.986001980626017 y.n=(y-a)*e^..= 4.119e-05 Inner loop: < 1 it --> adj=g*y= 4.119e-05 --Outer Loop: |tx - xinbta| > eps; tx= 0.999958501307107284 tx-xinbta = -4.119e-05 y=pbeta(x..)= 0.838616489670483 y.n=(y-a)*e^..= 0.0004112 Inner loop: < 1 it --> adj=g*y= 0.0004112 --Outer Loop: |tx - xinbta| > eps; tx= 0.999547280519047732 tx-xinbta = -0.0004112 y=pbeta(x..)= 0.501104895100778 y.n=(y-a)*e^..= 0.0009539 Inner loop: < 1 it --> adj=g*y= 0.0009539 --Outer Loop: |tx - xinbta| > eps; tx= 0.99859333619867674 tx-xinbta = -0.0009539 y=pbeta(x..)= 0.235563285738197 y.n=(y-a)*e^..= 0.001114 Inner loop: < 1 it --> adj=g*y= 0.001114 --Outer Loop: |tx - xinbta| > eps; tx= 0.997479445490663696 tx-xinbta = -0.001114 y=pbeta(x..)= 0.112234294496294 y.n=(y-a)*e^..= 0.0008728 Inner loop: < 1 it --> adj=g*y= 0.0008728 --Outer Loop: |tx - xinbta| > eps; tx= 0.996606669705101256 tx-xinbta = -0.0008728 y=pbeta(x..)= 0.0652991888096822 y.n=(y-a)*e^..= 0.0003853 Inner loop: < 1 it --> adj=g*y= 0.0003853 --Outer Loop: |tx - xinbta| > eps; tx= 0.996221360746289197 tx-xinbta = -0.0003853 y=pbeta(x..)= 0.0517477737754122 y.n=(y-a)*e^..= 5.634e-05 Inner loop: < 1 it --> adj=g*y= 5.634e-05 --Outer Loop: |tx - xinbta| > eps; tx= 0.99616502519828698 tx-xinbta = -5.634e-05 y=pbeta(x..)= 0.0500307755290968 y.n=(y-a)*e^..= 1.028e-06 Inner loop: < 1 it --> adj=g*y= 1.028e-06 --Outer Loop: |tx - xinbta| > eps; tx= 0.996163997249045807 tx-xinbta = -1.028e-06 y=pbeta(x..)= 0.0500000099834619 y.n=(y-a)*e^..= 3.337e-10 Inner loop: < 1 it --> adj=g*y= 3.337e-10 --Outer Loop: |tx - xinbta| > eps; tx= 0.996163996915366612 tx-xinbta = -3.337e-10 y=pbeta(x..)= 0.0500000000000011 y.n=(y-a)*e^..= 3.734e-17 Inner loop: < 1 it --> adj=g*y= 3.734e-17 10 outer iterations [1] 0.996164 > > showProc.time()#----------------------------------------------------------------------------- Time (user system elapsed): 0.04 0 0.05 > > DEBUG <- TRUE ## FIXME > > ## This is a "version" of qnormUappr() -- > ## but that has "regular" > ## > ## lower.tail=TRUE, log.p=FALSE > ## --- MM(2019) to MM(2008): I think this is a bad idea > ## ---> ../TODO: Use *argument* lower.tail=FALSE for all qnormU*() : "U" means "Upper tail" !!!! > qnormU <- function(u, lu) { + if(missing(u)) + qnorm(lu, lower.tail=FALSE,log.p=TRUE) + else + qnorm(u, lower.tail=FALSE) + } > > nln <- "\n----------------\n" > for(a in (10^(2:8))/2) + cat("p=", a, qbeta.R(.05, a, 1/2, low.bnd = 1e-9, qnormU = qnormU), + nln) p= 50 0.9621292 ---------------- p= 500 0.996164 ---------------- p= 5000 0.9996159 ---------------- p= 50000 0.9999616 ---------------- p= 5e+05 0.9999962 ---------------- p= 5e+06 0.9999996 ---------------- p= 5e+07 0.5 ---------------- > showProc.time() Time (user system elapsed): 0 0.01 0.02 > > > for(a in (10^(2:16))/2) cat("p=",a, qb <- qbeta.R(.05, a, 1/2, low.bnd = 0), 1-qb,nln) p= 50 0.9621292 0.03787078 ---------------- p= 500 0.996164 0.003836003 ---------------- p= 5000 0.9996159 0.0003840913 ---------------- p= 50000 0.9999616 3.841404e-05 ---------------- p= 5e+05 0.9999962 3.841453e-06 ---------------- p= 5e+06 0.9999996 3.841458e-07 ---------------- p= 5e+07 1 3.841459e-08 ---------------- p= 5e+08 0.5 0.5 ---------------- p= 5e+09 0.5 0.5 ---------------- p= 5e+10 0.5 0.5 ---------------- p= 5e+11 0.5 0.5 ---------------- p= 5e+12 0.5 0.5 ---------------- p= 5e+13 0.5 0.5 ---------------- p= 5e+14 0.5 0.5 ---------------- p= 5e+15 0.5 0.5 ---------------- > ## ---- -- SENSATIONAL ! -------- > a <- 5e13; cat("p=",a, qb <- qbeta.R(.05, a, 1/2, low.bnd = 0), 1-qb,nln) #-- problem: 0 outer it. p= 5e+13 0.5 0.5 ---------------- > a <- 5e14; cat("p=",a, qb <- qbeta.R(.05, a, 1/2, low.bnd = 0), 1-qb,nln) #-- problem: 19 outer it p= 5e+14 0.5 0.5 ---------------- > a <- 5e15; cat("p=",a, qb <- qbeta.R(.05, a, 1/2, low.bnd = 0), 1-qb,nln) #-- the last one "working" p= 5e+15 0.5 0.5 ---------------- > a <- 5e16; cat("p=",a, qb <- qbeta.R(.05, a, 1/2, low.bnd = 0), 1-qb,nln) # still fails: init xinbta =1 p= 5e+16 0.5 0.5 ---------------- > > for(a in (10^(2:16))/2) cat("p=",a, qb <- qbeta.R(.001, a, 1/2, low.bnd = 0), 1-qb,nln) p= 50 0.8968977 0.1031023 ---------------- p= 500 0.9892255 0.01077451 ---------------- p= 5000 0.9989178 0.001082225 ---------------- p= 50000 0.9998917 0.0001082703 ---------------- p= 5e+05 0.9999892 1.082751e-05 ---------------- p= 5e+06 0.9999989 1.082756e-06 ---------------- p= 5e+07 0.9999999 1.082757e-07 ---------------- p= 5e+08 1 1.082757e-08 ---------------- p= 5e+09 1 1.082757e-09 ---------------- p= 5e+10 1 1.082756e-10 ---------------- p= 5e+11 1 1.082667e-11 ---------------- p= 5e+12 1 1.082689e-12 ---------------- p= 5e+13 1 1.076916e-13 ---------------- p= 5e+14 1 9.325873e-15 ---------------- p= 5e+15 0.5 0.5 ---------------- > ## ---- > > ## BOTH p & q are BIG: --- sometimes LONG loop ! > for(a in (10^(2:15))/2) cat("p=",a, qb <- qbeta.R(.05, a, a/2, low.bnd = 0), 1-qb,nln) p= 50 0.5751702 0.4248298 ---------------- p= 500 0.6381189 0.3618811 ---------------- p= 5000 0.6576885 0.3423115 ---------------- p= 50000 0.6638328 0.3361672 ---------------- p= 5e+05 0.6657711 0.3342289 ---------------- p= 5e+06 0.6663835 0.3336165 ---------------- p= 5e+07 0.6665771 0.3334229 ---------------- p= 5e+08 0.6666384 0.3333616 ---------------- p= 5e+09 0.6666577 0.3333423 ---------------- p= 5e+10 0.6666638 0.3333362 ---------------- p= 5e+11 0.6666658 0.3333342 ---------------- p= 5e+12 0.6666664 0.3333336 ---------------- p= 5e+13 0.6666666 0.3333334 ---------------- p= 5e+14 0.6666666 0.3333334 ---------------- > ## ---- !WOW! > for(a in (10^(2:15))/2) cat("p=",a, qb <- qbeta.R(.05, a, a/2, low.bnd = 0, qnormU = qnormU, + f.acu = function(...)1e-26),1-qb,nln) p= 50 0.5751702 0.4248298 ---------------- p= 500 0.6381189 0.3618811 ---------------- p= 5000 0.6576885 0.3423115 ---------------- p= 50000 0.6638328 0.3361672 ---------------- p= 5e+05 0.6657711 0.3342289 ---------------- p= 5e+06 0.6663835 0.3336165 ---------------- p= 5e+07 0.6665771 0.3334229 ---------------- p= 5e+08 0.6666384 0.3333616 ---------------- p= 5e+09 0.6666577 0.3333423 ---------------- p= 5e+10 0.6666638 0.3333362 ---------------- p= 5e+11 0.6666658 0.3333342 ---------------- p= 5e+12 0.6666664 0.3333336 ---------------- p= 5e+13 0.6666666 0.3333334 ---------------- p= 5e+14 0.6666666 0.3333334 ---------------- > > options(op)# revert to usual > showProc.time() Time (user system elapsed): 0.02 0 0.01 > > > ## When using a new qbeta(.) C-function with low.bnd=0.0, up.bnd = 1.0 > ##--- STILL FAILS at 10^10 ... why ? ... > op <- options(warn=1) # immediate {there are only 2, now} > for(a in (10^(2:16))/2) + cat("p=",a, qb <- qbeta(.05, a, 1/2), 1-qb, " relE{pb, .05): ", formatC(pbeta(qb, a, 1/2)/.05 - 1), + nln) p= 50 0.9621292 0.03787078 relE{pb, .05): 3.109e-15 ---------------- p= 500 0.996164 0.003836003 relE{pb, .05): 2.243e-14 ---------------- p= 5000 0.9996159 0.0003840913 relE{pb, .05): 2.198e-13 ---------------- p= 50000 0.9999616 3.841404e-05 relE{pb, .05): -7.92e-13 ---------------- p= 5e+05 0.9999962 3.841453e-06 relE{pb, .05): -3.053e-11 ---------------- p= 5e+06 0.9999996 3.841458e-07 relE{pb, .05): 6.881e-11 ---------------- p= 5e+07 1 3.841459e-08 relE{pb, .05): 2.337e-09 ---------------- p= 5e+08 1 3.841459e-09 relE{pb, .05): -1.475e-09 ---------------- p= 5e+09 1 3.841459e-10 relE{pb, .05): -3.296e-07 ---------------- p= 5e+10 1 3.84146e-11 relE{pb, .05): -9.914e-07 ---------------- p= 5e+11 1 3.841483e-12 relE{pb, .05): -1.423e-05 ---------------- p= 5e+12 1 3.841372e-13 relE{pb, .05): 5.198e-05 ---------------- p= 5e+13 1 3.841372e-14 relE{pb, .05): 5.198e-05 ---------------- Warning in qbeta(0.05, a, 1/2) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.0020314 is not accurate p= 5e+14 1 3.885781e-15 relE{pb, .05): -0.02607 ---------------- Warning in qbeta(0.05, a, 1/2) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.017999 is not accurate p= 5e+15 1 3.330669e-16 relE{pb, .05): 0.36 ---------------- > options(op) > > showProc.time() Time (user system elapsed): 0.02 0 0.02 > > > ## Had tests for qbeta() built on AS 109 --- but now R *does* use AS 109 ! > ## This will probably all be irrelevant now, see also > ## ==> ~/R/MM/NUMERICS/dpq-functions/beta-gamma-etc/qbeta109.R > > > > pr.val <- c(.2,.4, 10^(-1:-5)); pr.val <- sort(c(pr.val, 1/2, 1-pr.val)) > > p.seq <- sort(outer(c(1,3),10^(3:16))) > p.seq <- c(2^(-15:6),10^(2:18)) > p.seq <- c(1e5^c(-4:-1), 10^c(-2:2,10:18)) > > ## q = 1 ---------- know exact formula ... > for(x in pr.val) { + cat("-----------\nx=", format01prec(x, digits = 6),":\n~~~~~~~~~~\n") + for(p in p.seq) { + t0 <- proc.time()[1] + qb <- qbeta(x, p, 1) # takes forever + C1 <- (t1 <- proc.time()[1]) - t0 + pb <- pbeta(x^(1/p), p, 1) + C2 <- (t2 <- proc.time()[1]) - t1 + fcat( + ## NO Debug: + paste0(formatC(p,w = 9), "|",formatC(C1, digits = 2, wid = 5),"|"), + ## In any case: + format01prec(qb, digits = 10), + qb - x^(1/p), qb /( x^(1/p))-1, pb - x + ##, f.dig=3, f.wid=8 # << Debug + ) + } + } ----------- x= 1e-05 : ~~~~~~~~~~ 1e-20| 0| 0 0 NaN -1e-05 1e-15| 0| 0 0 NaN -1e-05 1e-10| 0| 0 0 NaN -1e-05 1e-05| 0| 0 0 NaN -1e-05 0.01| 0| 0 0 NaN -1e-05 0.1| 0| 1e-50 7.834e-65 7.772e-15 -8.47e-21 1| 0| 1e-05 0 0 -3.388e-21 10| 0| 0.316227766 0 0 2.541e-20 100| 0| 0.8912509381 0 0 3.219e-20 1e+10| 0| 1-1.1513e-09 0 0 4.006e-12 1e+11| 0| 1-1.1513e-10 0 0 1.517e-11 1e+12| 0| 1-1.1513e-11 0 0 2.372e-10 1e+13| 0| 1-1.1513e-12 0 0 -8.73e-10 1e+14| 0| 1-1.1513e-13 0 0 -8.73e-10 1e+15| 0| 1-1.1435e-14 1.11e-16 2.22e-16 -3.284e-07 1e+16| 0| 1-1.1102e-15 0 0 5.079e-06 1e+17| 0| 1-1.1102e-16 0 0 5.079e-06 1e+18| 0| 1-1.1102e-16 -1.11e-16 -1.11e-16 1 ----------- x= 0.0001 : ~~~~~~~~~~ 1e-20| 0| 0 0 NaN -0.0001 1e-15| 0| 0 0 NaN -0.0001 1e-10| 0| 0 0 NaN -0.0001 1e-05| 0| 0 0 NaN -0.0001 0.01| 0| 0 0 NaN -0.0001 0.1| 0| 1e-40 6.525e-55 6.439e-15 -6.776e-20 1| 0| 0.0001 0 0 8.132e-20 10| 0| 0.3981071706 -5.551e-17 -1.11e-16 1.084e-19 100| 0| 0.9120108394 0 0 2.982e-19 1e+10| 0| 1-9.2103e-10 0 0 9.949e-12 1e+11| 0| 1-9.2103e-11 0 0 -3.227e-10 1e+12| 0| 1-9.2103e-12 0 0 4.118e-09 1e+13| 0| 1-9.2104e-13 0 0 -6.984e-09 1e+14| 0| 1-9.2149e-14 0 0 -4.501e-07 1e+15| 0| 1-9.1038e-15 1.11e-16 2.22e-16 -4.501e-07 1e+16| 0| 1-8.8818e-16 0 0 3.89e-05 1e+17| 0| 1-1.1102e-16 0 0 -8.492e-05 1e+18| 0| 1-1.1102e-16 -1.11e-16 -1.11e-16 0.9999 ----------- x= 0.001 : ~~~~~~~~~~ 1e-20| 0| 0 0 NaN -0.001 1e-15| 0| 0 0 NaN -0.001 1e-10| 0| 0 0 NaN -0.001 1e-05| 0| 0 0 NaN -0.001 0.01| 0| 1e-300 2.172e-314 2.176e-14 -2.168e-19 0.1| 0| 1e-30 6.481e-45 6.439e-15 -6.505e-19 1| 0| 0.001 0 0 2.168e-19 10| 0| 0.5011872336 0 0 -1.301e-18 100| 0| 0.9332543008 0 0 -3.686e-18 1e+10| 0| 1-6.9078e-10 0 0 -2.021e-10 1e+11| 0| 1-6.9078e-11 0 0 3.131e-09 1e+12| 0| 1-6.9078e-12 0 0 -5.238e-08 1e+13| 0| 1-6.9078e-13 0 0 -5.238e-08 1e+14| 0| 1-6.9056e-14 0 0 2.17e-06 1e+15| 0| 1-6.8834e-15 0 0 2.467e-05 1e+16| 0| 1-6.6613e-16 0 0 0.0002794 1e+17| 0| 1-1.1102e-16 0 0 -0.0009849 1e+18| 0| 1-1.1102e-16 -1.11e-16 -1.11e-16 0.999 ----------- x= 0.01 : ~~~~~~~~~~ 1e-20| 0| 0 0 NaN -0.01 1e-15| 0| 0 0 NaN -0.01 1e-10| 0| 0 0 NaN -0.01 1e-05| 0| 0 0 NaN -0.01 0.01| 0| 1e-200 1.537e-214 1.532e-14 -1.735e-18 0.1| 0| 1e-20 4.213e-35 4.219e-15 -3.469e-18 1| 0| 0.01 0 0 3.469e-18 10| 0| 0.6309573445 0 0 5.204e-18 100| 0| 0.954992586 0 0 4.857e-17 1e+10| 0| 1-4.6052e-10 0 0 -5.043e-09 1e+11| 0| 1-4.6052e-11 0 0 -1.614e-08 1e+12| 0| 1-4.6052e-12 0 0 -3.492e-07 1e+13| 0| 1-4.6052e-13 0 0 -3.492e-07 1e+14| 0| 1-4.6074e-14 0 0 -2.253e-05 1e+15| 0| 1-4.5519e-15 0 0 0.000547 1e+16| 0| 1-4.4409e-16 0 0 0.001785 1e+17| 0| 1-1.1102e-16 -1.11e-16 -1.11e-16 0.99 1e+18| 0| 1-1.1102e-16 -1.11e-16 -1.11e-16 0.99 ----------- x= 0.1 : ~~~~~~~~~~ 1e-20| 0| 0 0 NaN -0.1 1e-15| 0| 0 0 NaN -0.1 1e-10| 0| 0 0 NaN -0.1 1e-05| 0| 0 0 NaN -0.1 0.01| 0| 1e-100 -1.396e-115 -1.443e-15 0 0.1| 0| 1e-10 2.585e-25 2.665e-15 -2.776e-17 1| 0| 0.1 0 0 1.388e-17 10| 0| 0.7943282347 0 0 4.163e-17 100| 0| 0.977237221 0 0 -9.714e-17 1e+10| 0| 1-2.3026e-10 0 0 3.032e-08 1e+11| 0| 1-2.3026e-11 0 0 4.744e-07 1e+12| 0| 1-2.3026e-12 0 0 -1.746e-06 1e+13| 0| 1-2.3026e-13 0 0 -1.746e-06 1e+14| 0| 1-2.2982e-14 0 0 0.0004433 1e+15| 0| 1-2.3315e-15 0 0 -0.002847 1e+16| 0| 1-2.2204e-16 0 0 0.008561 1e+17| 0| 1-1.1102e-16 -1.11e-16 -1.11e-16 0.9 1e+18| 0| 1-1.1102e-16 -1.11e-16 -1.11e-16 0.9 ----------- x= 0.2 : ~~~~~~~~~~ 1e-20| 0| 0 0 NaN -0.2 1e-15| 0| 0 0 NaN -0.2 1e-10| 0| 0 0 NaN -0.2 1e-05| 0| 0 0 NaN -0.2 0.01| 0| 1.2676506e-70 -1.609e-86 -1.11e-16 0 0.1| 0| 1.024e-07 2.647e-22 2.665e-15 -5.551e-17 1| 0| 0.2 0 0 0 10| 0| 0.8513399225 0 0 -8.327e-17 100| 0| 0.9840344434 0 0 -5.274e-16 1e+10| 0| 1-1.6094e-10 0 0 -4.52e-08 1e+11| 0| 1-1.6094e-11 0 0 6.21e-07 1e+12| 0| 1-1.6095e-12 0 0 -1.048e-05 1e+13| 0| 1-1.6098e-13 0 0 -7.708e-05 1e+14| 0| 1-1.6098e-14 0 0 -7.708e-05 1e+15| 0| 1-1.5543e-15 0 0 0.01133 1e+16| 0| 1-1.1102e-16 0 0 0.1295 1e+17| 0| 1-1.1102e-16 -1.11e-16 -1.11e-16 0.8 1e+18| 0| 1-1.1102e-16 -1.11e-16 -1.11e-16 0.8 ----------- x= 0.4 : ~~~~~~~~~~ 1e-20| 0| 0 0 NaN -0.4 1e-15| 0| 0 0 NaN -0.4 1e-10| 0| 0 0 NaN -0.4 1e-05| 0| 0 0 NaN -0.4 0.01| 0| 1.606938044e-40 1.427e-55 8.882e-16 0 0.1| 0| 0.0001048576 2.304e-19 2.22e-15 -5.551e-17 1| 0| 0.4 0 0 0 10| 0| 0.9124435366 0 0 -2.22e-16 100| 0| 0.9908789442 0 0 -5.551e-16 1e+10| 0| 1-9.1629e-11 0 0 1.42e-07 1e+11| 0| 1-9.1629e-12 0 0 5.861e-07 1e+12| 0| 1-9.1627e-13 0 0 9.468e-06 1e+13| 0| 1-9.1593e-14 0 0 0.0001427 1e+14| 0| 1-9.2149e-15 0 0 -0.002072 1e+15| 0| 1-8.8818e-16 0 0 0.0114 1e+16| 0| 1-1.1102e-16 0 0 -0.07051 1e+17| 0| 1-1.1102e-16 -1.11e-16 -1.11e-16 0.6 1e+18| 0| 1-1.1102e-16 -1.11e-16 -1.11e-16 0.6 ----------- x= 0.5 : ~~~~~~~~~~ 1e-20| 0| 0 0 NaN -0.5 1e-15| 0| 5.562684646e-309 5.563e-309 Inf -0.5 1e-10| 0| 5.562684646e-309 5.563e-309 Inf -0.5 1e-05| 0| 5.562684646e-309 5.563e-309 Inf -0.5 0.01| 0| 7.888609052e-31 -1.314e-45 -1.665e-15 0 0.1| 0| 0.0009765625 1.952e-18 1.998e-15 -5.551e-17 1| 0| 0.5 0 0 0 10| 0| 0.9330329915 0 0 -1.11e-16 100| 0| 0.9930924954 0 0 8.882e-16 1e+10| 0| 1-6.9315e-11 0 0 2.647e-07 1e+11| 0| 1-6.9315e-12 0 0 8.198e-07 1e+12| 0| 1-6.9311e-13 0 0 1.747e-05 1e+13| 0| 1-6.9278e-14 0 0 0.000184 1e+14| 0| 1-6.8834e-15 0 0 0.00241 1e+15| 0| 1-6.6613e-16 0 0 0.01369 1e+16| 0| 1-1.1102e-16 0 0 -0.1705 1e+17| 0| 1-1.1102e-16 -1.11e-16 -1.11e-16 0.5 1e+18| 0| 1-1.1102e-16 -1.11e-16 -1.11e-16 0.5 ----------- x= 0.6 : ~~~~~~~~~~ 1e-20| 0| 5.562684646e-309 5.563e-309 Inf -0.6 1e-15| 0| 5.562684646e-309 5.563e-309 Inf -0.6 1e-10| 0| 5.562684646e-309 5.563e-309 Inf -0.6 1e-05| 0| 5.562684646e-309 5.563e-309 Inf -0.6 0.01| 0| 6.533186235e-23 4.702e-38 6.661e-16 0 0.1| 0| 0.0060466176 8.674e-18 1.332e-15 -1.11e-16 1| 0| 0.6 0 0 0 10| 0| 0.9502002165 0 0 2.22e-16 100| 0| 0.9949047687 0 0 -2.998e-15 1e+10| 0| 1-5.1083e-11 0 0 -1.214e-07 1e+11| 0| 1-5.1082e-12 0 0 5.447e-07 1e+12| 0| 1-5.1081e-13 0 0 7.206e-06 1e+13| 0| 1-5.107e-14 0 0 7.382e-05 1e+14| 0| 1-5.107e-15 0 0 7.382e-05 1e+15| 0| 1-5.5511e-16 0 0 -0.02599 1e+16| 0| 1-0 0 0 0.4 1e+17| 0| 1-0 0 0 0.4 1e+18| 0| 1-0 0 0 0.4 ----------- x= 0.8 : ~~~~~~~~~~ 1e-20| 0| 5.562684646e-309 5.563e-309 Inf -0.8 1e-15| 0| 5.562684646e-309 5.563e-309 Inf -0.8 1e-10| 0| 5.562684646e-309 5.563e-309 Inf -0.8 1e-05| 0| 5.562684646e-309 5.563e-309 Inf -0.8 0.01| 0| 2.037035976e-10 -2.068e-25 -9.992e-16 0 0.1| 0| 0.1073741824 1.249e-16 1.11e-15 -1.11e-16 1| 0| 0.8 0 0 0 10| 0| 0.9779327685 0 0 2.22e-16 100| 0| 0.9977710523 0 0 -2.331e-15 1e+10| 0| 1-2.2314e-11 0 0 -1.395e-07 1e+11| 0| 1-2.2314e-12 0 0 -1.395e-07 1e+12| 0| 1-2.2315e-13 0 0 -9.021e-06 1e+13| 0| 1-2.2315e-14 0 0 -9.021e-06 1e+14| 0| 1-2.2204e-15 0 0 0.0008796 1e+15| 0| 1-2.2204e-16 0 0 0.0008796 1e+16| 0| 1-0 0 0 0.2 1e+17| 0| 1-0 0 0 0.2 1e+18| 0| 1-0 0 0 0.2 ----------- x= 0.9 : ~~~~~~~~~~ 1e-20| 0| 5.562684646e-309 5.563e-309 Inf -0.9 1e-15| 0| 5.562684646e-309 5.563e-309 Inf -0.9 1e-10| 0| 5.562684646e-309 5.563e-309 Inf -0.9 1e-05| 0| 5.562684646e-309 5.563e-309 Inf -0.9 0.01| 0| 2.656139889e-05 2.711e-20 1.11e-15 0 0.1| 0| 0.3486784401 5.551e-17 2.22e-16 -1.11e-16 1| 0| 0.9 0 0 0 10| 0| 0.9895192582 0 0 4.441e-16 100| 0| 0.9989469497 0 0 5.551e-16 1e+10| 0| 1-1.0536e-11 0 0 3.156e-07 1e+11| 0| 1-1.0536e-12 0 0 3.156e-07 1e+12| 0| 1-1.0536e-13 0 0 3.156e-07 1e+13| 0| 1-1.0547e-14 0 0 -9.96e-05 1e+14| 0| 1-9.992e-16 0 0 0.00491 1e+15| 0| 1-1.1102e-16 0 0 -0.005081 1e+16| 0| 1-0 0 0 0.1 1e+17| 0| 1-0 0 0 0.1 1e+18| 0| 1-0 0 0 0.1 ----------- x= 0.99 : ~~~~~~~~~~ 1e-20| 0| 5.562684646e-309 5.563e-309 Inf -0.99 1e-15| 0| 5.562684646e-309 5.563e-309 Inf -0.99 1e-10| 0| 5.562684646e-309 5.563e-309 Inf -0.99 1e-05| 0| 5.562684646e-309 5.563e-309 Inf -0.99 0.01| 0| 0.3660323413 -1.665e-16 -4.441e-16 0 0.1| 0| 0.904382075 0 0 0 1| 0| 0.99 0 0 0 10| 0| 0.9989954713 0 0 -1.11e-16 100| 0| 0.9998995017 0 0 8.882e-16 1e+10| 0| 1-1.0051e-12 0 0 -5.081e-07 1e+11| 0| 1-1.0048e-13 0 0 2.789e-06 1e+12| 0| 1-1.0103e-14 0 0 -5.217e-05 1e+13| 0| 1-9.992e-16 0 0 5.775e-05 1e+14| 0| 1-1.1102e-16 0 0 -0.001041 1e+15| 0| 1-0 0 0 0.01 1e+16| 0| 1-0 0 0 0.01 1e+17| 0| 1-0 0 0 0.01 1e+18| 0| 1-0 0 0 0.01 ----------- x= 0.999 : ~~~~~~~~~~ 1e-20| 0| 5.562684646e-309 5.563e-309 Inf -0.999 1e-15| 0| 5.562684646e-309 5.563e-309 Inf -0.999 1e-10| 0| 5.562684646e-309 5.563e-309 Inf -0.999 1e-05| 0| 3.538527688e-44 -2.763e-57 -7.805e-14 0 0.01| 0| 0.9047921471 0 0 0 0.1| 0| 0.9900448802 0 0 0 1| 0| 0.999 0 0 0 10| 0| 0.999899955 0 0 0 100| 0| 0.999989995 0 0 -1.11e-16 1e+10| 0| 1-1.0003e-13 0 0 1.892e-07 1e+11| 0| 1-9.992e-15 0 0 1.298e-06 1e+12| 0| 1-9.992e-16 0 0 1.298e-06 1e+13| 0| 1-1.1102e-16 0 0 -0.0001096 1e+14| 0| 1-0 0 0 0.001 1e+15| 0| 1-0 0 0 0.001 1e+16| 0| 1-0 0 0 0.001 1e+17| 0| 1-0 0 0 0.001 1e+18| 0| 1-0 0 0 0.001 ----------- x= 0.9999 : ~~~~~~~~~~ 1e-20| 0| 5.562684646e-309 5.563e-309 Inf -0.9999 1e-15| 0| 5.562684646e-309 5.563e-309 Inf -0.9999 1e-10| 0| 5.562684646e-309 5.563e-309 Inf -0.9999 1e-05| 0| 4.537723396e-05 -2.643e-19 -5.773e-15 0 0.01| 0| 0.9900493387 0 0 0 0.1| 0| 0.9990004499 0 0 0 1| 0| 0.9999 0 0 0 10| 0| 0.9999899995 0 0 -2.22e-16 100| 0| 0.999999 0 0 -3.553e-15 1e+10| 0| 1-9.992e-15 0 0 8.492e-08 1e+11| 0| 1-9.992e-16 0 0 8.492e-08 1e+12| 0| 1-1.1102e-16 0 0 -1.102e-05 1e+13| 0| 1-0 0 0 1e-04 1e+14| 0| 1-0 0 0 1e-04 1e+15| 0| 1-0 0 0 1e-04 1e+16| 0| 1-0 0 0 1e-04 1e+17| 0| 1-0 0 0 1e-04 1e+18| 0| 1-0 0 0 1e-04 ----------- x= 0.99999 : ~~~~~~~~~~ 1e-20| 0| 5.562684646e-309 5.563e-309 Inf -1 1e-15| 0| 5.562684646e-309 5.563e-309 Inf -1 1e-10| 0| 5.562684646e-309 5.563e-309 Inf -1 1e-05| 0| 0.3678776018 -2.22e-16 -5.551e-16 0 0.01| 0| 0.9990004948 0 0 0 0.1| 0| 0.9999000045 0 0 0 1| 0| 0.99999 0 0 0 10| 0| 0.999999 0 0 3.331e-16 100| 0| 1-1e-07 0 0 7.772e-16 1e+10| 0| 1-9.992e-16 0 0 8.043e-09 1e+11| 0| 1-1.1102e-16 0 0 -1.102e-06 1e+12| 0| 1-0 0 0 1e-05 1e+13| 0| 1-0 0 0 1e-05 1e+14| 0| 1-0 0 0 1e-05 1e+15| 0| 1-0 0 0 1e-05 1e+16| 0| 1-0 0 0 1e-05 1e+17| 0| 1-0 0 0 1e-05 1e+18| 0| 1-0 0 0 1e-05 There were 31 warnings (use warnings() to see them) > > summary(warnings()) Summary of (a total of 31) warning messages: 1x : In qbeta(x, p, 1) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.00017936 is not accurate 1x : In qbeta(x, p, 1) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.00027943 is not accurate 2x : In qbeta(x, p, 1) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.001 is not accurate 1x : In qbeta(x, p, 1) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.00013269 is not accurate 1x : In qbeta(x, p, 1) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.0005613 is not accurate 1x : In qbeta(x, p, 1) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.0017854 is not accurate 2x : In qbeta(x, p, 1) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.01 is not accurate 1x : In qbeta(x, p, 1) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.0026986 is not accurate 1x : In qbeta(x, p, 1) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.002847 is not accurate 1x : In qbeta(x, p, 1) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.0085607 is not accurate 2x : In qbeta(x, p, 1) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.1 is not accurate 1x : In qbeta(x, p, 1) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.0021549 is not accurate 1x : In qbeta(x, p, 1) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.011335 is not accurate 1x : In qbeta(x, p, 1) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.091439 is not accurate 2x : In qbeta(x, p, 1) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.2 is not accurate 1x : In qbeta(x, p, 1) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.0016303 is not accurate 1x : In qbeta(x, p, 1) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.0020724 is not accurate 1x : In qbeta(x, p, 1) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.011404 is not accurate 1x : In qbeta(x, p, 1) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.29144 is not accurate 2x : In qbeta(x, p, 1) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.4 is not accurate 1x : In qbeta(x, p, 1) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.0047889 is not accurate 1x : In qbeta(x, p, 1) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.0031368 is not accurate 1x : In qbeta(x, p, 1) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.013691 is not accurate 1x : In qbeta(x, p, 1) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.39144 is not accurate 2x : In qbeta(x, p, 1) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.5 is not accurate > showProc.time() Time (user system elapsed): 0.1 0 0.11 > > ## q = 1/2 (T-distrib.) > ## qt is defined as > ## qt(x,n) := - sqrt((1 / qbeta( 2*x , n/2, 1/2) - 1) * n) ## for x <= 1/2 > ## qt(x,n) := sqrt((1 / qbeta( 2*(1-x), n/2, 1/2) - 1) * n) ## for x > 1/2 > for(x in pr.val) { + cat("-----------\nx=", format01prec(x, digits = 6),":\n~~~~~~~~~~\n") + for(p in p.seq[p.seq >= 1/2]) { + qb0 <- 1/(1+ qt(x/2, df = 2*p)^2/(2*p)) + t0 <- proc.time()[1] + qb <- qbeta(x, p, 1/2) # takes forever + C1 <- (t1 <- proc.time()[1]) - t0 + pb <- pbeta(qb, p, 1/2) + C2 <- (t2 <- proc.time()[1]) - t1 + ## NO Debug: + fcat(p,paste("|",t2-t0,"|",sep = ''),format01prec(c(qb,qb0),digits = 10), + ## Debug: fcat(format01prec(qb,digits=10), + paste("|",formatC(qb0-qb,dig = 8,w = 12)), + qb/qb0 - 1, pb/x - 1, f.wid = 12) + } + } ----------- x= 1e-05 : ~~~~~~~~~~ 1 |0| 1.99999e-05 1.99999e-05 | -6.7762636e-21 4.441e-16 0 10 |0| 0.3685587802 0.3685587802 | -5.5511151e-17 2.22e-16 -2.554e-15 100 |0| 0.9068308826 0.9068308826 | 1.110223e-16 -1.11e-16 6.661e-16 1e+10 |0| 1-9.7557e-10 1-9.7557e-10 | 1.110223e-16 -1.11e-16 -2.156e-07 1e+11 |0| 1-9.7557e-11 1-9.7557e-11 | 0 0 3.276e-06 1e+12 |0| 1-9.7558e-12 1-9.7558e-12 | 0 0 -4.322e-05 1e+13 |0| 1-9.7555e-13 1-9.7566e-13 | -1.110223e-16 2.22e-16 0.0001893 1e+14 |0| 1-9.7589e-14 1-9.7478e-14 | 1.110223e-16 -1.11e-16 -0.003292 1e+15 |0| 1-9.4369e-15 1-9.77e-15 | -3.3306691e-16 4.441e-16 0.3966 1e+16 |0| 1-7.7716e-16 1-8.8818e-16 | -1.110223e-16 2.22e-16 7.064 1e+17 |0| 1-1.1102e-16 1-0 | 1.110223e-16 -1.11e-16 -0.7549 1e+18 |0| 1-1.1102e-16 1-0 | 1.110223e-16 -1.11e-16 -1 ----------- x= 0.0001 : ~~~~~~~~~~ 1 |0| 0.00019999 0.00019999 | 2.7105054e-20 -1.11e-16 0 10 |0| 0.4608349887 0.4608349887 | 0 0 4.441e-16 100 |0| 0.926935061 0.926935061 | -1.110223e-16 2.22e-16 1.554e-15 1e+10 |0| 1-7.5684e-10 1-7.5684e-10 | 0 0 -6.885e-10 1e+11 |0| 1-7.5684e-11 1-7.5683e-11 | 1.110223e-16 -1.11e-16 -4.702e-06 1e+12 |0| 1-7.5684e-12 1-7.5684e-12 | 0 0 -3.998e-05 1e+13 |0| 1-7.5684e-13 1-7.5673e-13 | 1.110223e-16 -1.11e-16 -3.998e-05 1e+14 |0| 1-7.5717e-14 1-7.5717e-14 | 0 0 -0.003562 1e+15 |0| 1-7.3275e-15 1-7.5495e-15 | -2.220446e-16 2.22e-16 0.291 1e+16 |0| 1-7.7716e-16 1-6.6613e-16 | 1.110223e-16 -1.11e-16 -0.1936 1e+17 |0| 1-1.1102e-16 1-0 | 1.110223e-16 -1.11e-16 -0.9755 1e+18 |0| 1-1.1102e-16 1-0 | 1.110223e-16 -1.11e-16 -1 ----------- x= 0.001 : ~~~~~~~~~~ 1 |0| 0.001999 0.001999 | 4.3368087e-19 -2.22e-16 0 10 |0| 0.5744027377 0.5744027377 | 0 0 -8.882e-16 100 |0| 0.9471737495 0.9471737495 | 0 0 1.11e-15 1e+10 |0| 1-5.4138e-10 1-5.4138e-10 | 0 0 3.549e-07 1e+11 |0| 1-5.4138e-11 1-5.4138e-11 | 0 0 2.755e-06 1e+12 |0| 1-5.4138e-12 1-5.4139e-12 | -1.110223e-16 2.22e-16 2.755e-06 1e+13 |0.02| 1-5.4134e-13 1-5.4134e-13 | 0 0 0.0003626 1e+14 |0| 1-5.4179e-14 1-5.4179e-14 | 0 0 -0.004425 1e+15 |0| 1-5.4401e-15 1-5.3291e-15 | 1.110223e-16 -1.11e-16 -0.02801 1e+16 |0| 1-5.5511e-16 1-4.4409e-16 | 1.110223e-16 -1.11e-16 -0.1378 1e+17 |0| 1-1.1102e-16 1-0 | 1.110223e-16 -1.11e-16 -0.9975 1e+18 |0| 1-1.1102e-16 1-0 | 1.110223e-16 -1.11e-16 -1 ----------- x= 0.01 : ~~~~~~~~~~ 1 |0| 0.0199 0.0199 | 0 0 0 10 |0| 0.711846165 0.711846165 | 1.110223e-16 -1.11e-16 -3.331e-16 100 |0| 0.967289652 0.967289652 | 2.220446e-16 -2.22e-16 5.329e-15 1e+10 |0| 1-3.3174e-10 1-3.3174e-10 | 0 0 -2.673e-07 1e+11 |0| 1-3.3174e-11 1-3.3175e-11 | -1.110223e-16 2.22e-16 2.226e-06 1e+12 |0| 1-3.3175e-12 1-3.3173e-12 | 1.110223e-16 -1.11e-16 -1.024e-05 1e+13 |0| 1-3.3173e-13 1-3.3173e-13 | 0 0 0.0001144 1e+14 |0| 1-3.3196e-14 1-3.3085e-14 | 1.110223e-16 -1.11e-16 -0.002376 1e+15 |0| 1-3.3307e-15 1-3.3307e-15 | 0 0 -0.01473 1e+16 |0| 1-3.3307e-16 1-2.2204e-16 | 1.110223e-16 -1.11e-16 -0.01473 1e+17 |0| 1-1.1102e-16 1-0 | 1.110223e-16 -1.11e-16 -0.9998 1e+18 |0| 1-1.1102e-16 1-0 | 1.110223e-16 -1.11e-16 -1 ----------- x= 0.1 : ~~~~~~~~~~ 1 |0| 0.19 0.19 | -8.3266727e-17 4.441e-16 0 10 |0| 0.8705245727 0.8705245727 | 0 0 8.882e-16 100 |0| 0.9865300031 0.9865300031 | 1.110223e-16 -1.11e-16 -6.217e-15 1e+10 |0| 1-1.3528e-10 1-1.3528e-10 | 0 0 6.253e-07 1e+11 |0| 1-1.3528e-11 1-1.3528e-11 | 1.110223e-16 -1.11e-16 -2.159e-06 1e+12 |0| 1-1.3528e-12 1-1.3527e-12 | 1.110223e-16 -1.11e-16 -4.393e-05 1e+13 |0| 1-1.3523e-13 1-1.3523e-13 | 0 0 0.0006524 1e+14 |0| 1-1.3545e-14 1-1.3545e-14 | 0 0 -0.00213 1e+15 |0| 1-1.3323e-15 1-1.3323e-15 | 0 0 0.02608 1e+16 |0| 1-1.1102e-16 1-2.2204e-16 | -1.110223e-16 2.22e-16 0.3619 1e+17 |0| 1-1.1102e-16 1-0 | 1.110223e-16 -1.11e-16 -1 1e+18 |0| 1-1.1102e-16 1-0 | 1.110223e-16 -1.11e-16 -1 ----------- x= 0.2 : ~~~~~~~~~~ 1 |0| 0.36 0.36 | 0 0 0 10 |0| 0.9192643242 0.9192643242 | 0 0 2.22e-16 100 |0| 0.9918013803 0.9918013803 | 0 0 -6.106e-15 1e+10 |0| 1-8.2119e-11 1-8.2119e-11 | 0 0 -4.859e-07 1e+11 |0| 1-8.2119e-12 1-8.2119e-12 | 0 0 -4.859e-07 1e+12 |0| 1-8.2123e-13 1-8.2112e-13 | 1.110223e-16 -1.11e-16 -6.13e-05 1e+13 |0| 1-8.2157e-14 1-8.2157e-14 | 0 0 -0.0005173 1e+14 |0| 1-8.2157e-15 1-8.2157e-15 | 0 0 -0.0005173 1e+15 |0| 1-7.7716e-16 1-8.8818e-16 | -1.110223e-16 2.22e-16 0.0625 1e+16 |0| 1-1.1102e-16 1-0 | 1.110223e-16 -1.11e-16 -0.319 1e+17 |0| 1-1.1102e-16 1-0 | 1.110223e-16 -1.11e-16 -1 1e+18 |0| 1-1.1102e-16 1-0 | 1.110223e-16 -1.11e-16 -1 ----------- x= 0.4 : ~~~~~~~~~~ 1 |0| 0.64 0.64 | 2.220446e-16 -3.331e-16 0 10 |0| 0.9643415968 0.9643415968 | 1.110223e-16 -1.11e-16 -2.22e-16 100 |0| 0.9964558016 0.9964558016 | 1.110223e-16 -1.11e-16 -6.106e-15 1e+10 |0| 1-3.5416e-11 1-3.5416e-11 | 0 0 -3.574e-07 1e+11 |0| 1-3.5416e-12 1-3.5416e-12 | 0 0 3.336e-06 1e+12 |0| 1-3.5416e-13 1-3.5416e-13 | 0 0 3.336e-06 1e+13 |0| 1-3.5416e-14 1-3.5527e-14 | -1.110223e-16 2.22e-16 3.336e-06 1e+14 |0| 1-3.5527e-15 1-3.5527e-15 | 0 0 -0.001841 1e+15 |0| 1-3.3307e-16 1-4.4409e-16 | -1.110223e-16 2.22e-16 0.03601 1e+16 |0| 1-1.1102e-16 1-0 | 1.110223e-16 -1.11e-16 -0.6595 1e+17 |0| 1-1.1102e-16 1-0 | 1.110223e-16 -1.11e-16 -1 1e+18 |0| 1-1.1102e-16 1-0 | 1.110223e-16 -1.11e-16 -1 ----------- x= 0.5 : ~~~~~~~~~~ 1 |0| 0.75 0.75 | 0 0 0 10 |0| 0.9769485817 0.9769485817 | 0 0 4.441e-16 100 |0| 0.9977222233 0.9977222233 | 0 0 -6.661e-16 1e+10 |0| 1-2.2747e-11 1-2.2747e-11 | -1.110223e-16 2.22e-16 3.15e-07 1e+11 |0| 1-2.2747e-12 1-2.2746e-12 | 1.110223e-16 -1.11e-16 -1.015e-05 1e+12 |0| 1-2.2748e-13 1-2.2737e-13 | 1.110223e-16 -1.11e-16 -3.107e-05 1e+13 |0| 1-2.276e-14 1-2.2649e-14 | 1.110223e-16 -1.11e-16 -0.0002402 1e+14 |0| 1-2.2204e-15 1-2.2204e-15 | 0 0 0.01031 1e+15 |0| 1-2.2204e-16 1-2.2204e-16 | 0 0 0.01031 1e+16 |0| 1-1.1102e-16 1-0 | 1.110223e-16 -1.11e-16 -0.7276 1e+17 |0| 1-1.1102e-16 1-0 | 1.110223e-16 -1.11e-16 -1 1e+18 |0| 1-1.1102e-16 1-0 | 1.110223e-16 -1.11e-16 -1 ----------- x= 0.6 : ~~~~~~~~~~ 1 |0| 0.84 0.84 | 0 0 0 10 |0| 0.9860015995 0.9860015995 | -1.110223e-16 2.22e-16 -5.551e-16 100 |0| 0.9986225288 0.9986225288 | 1.110223e-16 -1.11e-16 -8.105e-15 1e+10 |0| 1-1.375e-11 1-1.375e-11 | -1.110223e-16 2.22e-16 3.491e-07 1e+11 |0| 1-1.375e-12 1-1.3749e-12 | 1.110223e-16 -1.11e-16 -7.012e-06 1e+12 |0| 1-1.3745e-13 1-1.3745e-13 | 0 0 0.0001157 1e+13 |0| 1-1.3767e-14 1-1.3767e-14 | 0 0 -0.0003749 1e+14 |0| 1-1.3323e-15 1-1.3323e-15 | 0 0 0.009535 1e+15 |0| 1-1.1102e-16 1-2.2204e-16 | -1.110223e-16 2.22e-16 0.06248 1e+16 |0| 1-0 1-0 | 0 0 0.6667 1e+17 |0| 1-0 1-0 | 0 0 0.6667 1e+18 |0| 1-0 1-0 | 0 0 0.6667 ----------- x= 0.8 : ~~~~~~~~~~ 1 |0| 0.96 0.96 | 0 0 -1.11e-16 10 |0| 0.9967149848 0.9967149848 | 0 0 -9.992e-16 100 |0| 0.9996783247 0.9996783247 | 0 0 -1.632e-14 1e+10 |0| 1-3.2092e-12 1-3.2092e-12 | 0 0 1.032e-06 1e+11 |0| 1-3.2097e-13 1-3.2085e-13 | 1.110223e-16 -1.11e-16 -1.59e-05 1e+12 |0| 1-3.2085e-14 1-3.2196e-14 | -1.110223e-16 2.22e-16 2.643e-05 1e+13 |0| 1-3.2196e-15 1-3.1086e-15 | 1.110223e-16 -1.11e-16 -0.0003965 1e+14 |0| 1-3.3307e-16 1-2.2204e-16 | 1.110223e-16 -1.11e-16 -0.004584 1e+15 |0| 1-0 1-0 | 0 0 0.25 1e+16 |0| 1-0 1-0 | 0 0 0.25 1e+17 |0| 1-0 1-0 | 0 0 0.25 1e+18 |0| 1-0 1-0 | 0 0 0.25 ----------- x= 0.9 : ~~~~~~~~~~ 1 |0| 0.99 0.99 | 1.110223e-16 -1.11e-16 -1.11e-16 10 |0| 0.9991908114 0.9991908114 | 0 0 -6.661e-16 100 |0| 0.9999208516 0.9999208516 | 0 0 1.865e-14 1e+10 |0| 1-7.8959e-13 1-7.8959e-13 | 0 0 -3.633e-06 1e+11 |0| 1-7.8937e-14 1-7.9048e-14 | -1.110223e-16 2.22e-16 1.191e-05 1e+12 |0| 1-7.8826e-15 1-7.9936e-15 | -1.110223e-16 2.22e-16 8.966e-05 1e+13 |0| 1-7.7716e-16 1-8.8818e-16 | -1.110223e-16 2.22e-16 0.0008702 1e+14 |0| 1-1.1102e-16 1-0 | 1.110223e-16 -1.11e-16 -0.02051 1e+15 |0| 1-0 1-0 | 0 0 0.1111 1e+16 |0| 1-0 1-0 | 0 0 0.1111 1e+17 |0| 1-0 1-0 | 0 0 0.1111 1e+18 |0| 1-0 1-0 | 0 0 0.1111 ----------- x= 0.99 : ~~~~~~~~~~ 1 |0| 0.9999 0.9999 | 1.110223e-16 -1.11e-16 6.661e-16 10 |0| 0.9999919469 0.9999919469 | 0 0 -3.164e-14 100 |0| 1-7.8741e-07 1-7.8741e-07 | 0 0 1.454e-13 1e+10 |0| 1-7.8826e-15 1-7.7716e-15 | 1.110223e-16 -1.11e-16 -1.811e-05 1e+11 |0| 1-7.7716e-16 1-8.8818e-16 | -1.110223e-16 2.22e-16 5.34e-05 1e+12 |0| 1-1.1102e-16 1-0 | 1.110223e-16 -1.11e-16 -0.001908 1e+13 |0| 1-0 1-0 | 0 0 0.0101 1e+14 |0| 1-0 1-0 | 0 0 0.0101 1e+15 |0| 1-0 1-0 | 0 0 0.0101 1e+16 |0| 1-0 1-0 | 0 0 0.0101 1e+17 |0| 1-0 1-0 | 0 0 0.0101 1e+18 |0| 1-0 1-0 | 0 0 0.0101 ----------- x= 0.999 : ~~~~~~~~~~ 1 |0| 0.999999 0.999999 | 0 0 -1.432e-14 10 |0| 1-8.0527e-08 1-8.0527e-08 | -1.110223e-16 2.22e-16 -1.545e-13 100 |0| 1-7.8736e-09 1-7.8736e-09 | 0 0 -1.174e-12 1e+10 |0| 1-1.1102e-16 1-0 | 1.110223e-16 -1.11e-16 -0.0001891 1e+11 |0| 1-0 1-0 | 0 0 0.001001 1e+12 |0| 1-0 1-0 | 0 0 0.001001 1e+13 |0| 1-0 1-0 | 0 0 0.001001 1e+14 |0| 1-0 1-0 | 0 0 0.001001 1e+15 |0| 1-0 1-0 | 0 0 0.001001 1e+16 |0| 1-0 1-0 | 0 0 0.001001 1e+17 |0| 1-0 1-0 | 0 0 0.001001 1e+18 |0| 1-0 1-0 | 0 0 0.001001 ----------- x= 0.9999 : ~~~~~~~~~~ 1 |0| 1-1e-08 1-1e-08 | 0 0 -2.512e-13 10 |0| 1-8.0527e-10 1-8.0527e-10 | 0 0 -2.494e-12 100 |0| 1-7.8736e-11 1-7.8736e-11 | 1.110223e-16 -1.11e-16 -3.219e-11 1e+10 |0| 1-0 1-0 | 0 0 0.0001 1e+11 |0| 1-0 1-0 | 0 0 0.0001 1e+12 |0| 1-0 1-0 | 0 0 0.0001 1e+13 |0| 1-0 1-0 | 0 0 0.0001 1e+14 |0| 1-0 1-0 | 0 0 0.0001 1e+15 |0| 1-0 1-0 | 0 0 0.0001 1e+16 |0| 1-0 1-0 | 0 0 0.0001 1e+17 |0| 1-0 1-0 | 0 0 0.0001 1e+18 |0| 1-0 1-0 | 0 0 0.0001 ----------- x= 0.99999 : ~~~~~~~~~~ 1 |0| 1-1e-10 1-1e-10 | 0 0 -4.138e-13 10 |0| 1-8.0527e-12 1-8.0527e-12 | 0 0 3.282e-11 100 |0| 1-7.8737e-13 1-7.8737e-13 | 0 0 -3.85e-11 1e+10 |0| 1-0 1-0 | 0 0 1e-05 1e+11 |0| 1-0 1-0 | 0 0 1e-05 1e+12 |0| 1-0 1-0 | 0 0 1e-05 1e+13 |0| 1-0 1-0 | 0 0 1e-05 1e+14 |0| 1-0 1-0 | 0 0 1e-05 1e+15 |0| 1-0 1-0 | 0 0 1e-05 1e+16 |0| 1-0 1-0 | 0 0 1e-05 1e+17 |0| 1-0 1-0 | 0 0 1e-05 1e+18 |0| 1-0 1-0 | 0 0 1e-05 There were 32 warnings (use warnings() to see them) > > summary(warnings()) Summary of (a total of 32) warning messages: 2x : In qbeta(x, p, 1/2) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.00013776 is not accurate 2x : In qbeta(x, p, 1/2) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.001 is not accurate 1x : In qbeta(x, p, 1/2) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.00022814 is not accurate 2x : In qbeta(x, p, 1/2) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.00014731 is not accurate 2x : In qbeta(x, p, 1/2) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.01 is not accurate 1x : In qbeta(x, p, 1/2) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.0040508 is not accurate 1x : In qbeta(x, p, 1/2) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.0026078 is not accurate 1x : In qbeta(x, p, 1/2) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.064912 is not accurate 2x : In qbeta(x, p, 1/2) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.1 is not accurate 1x : In qbeta(x, p, 1/2) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.00010345 is not accurate 1x : In qbeta(x, p, 1/2) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.0029628 is not accurate 1x : In qbeta(x, p, 1/2) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.0125 is not accurate 1x : In qbeta(x, p, 1/2) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.16491 is not accurate 2x : In qbeta(x, p, 1/2) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.2 is not accurate 1x : In qbeta(x, p, 1/2) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.0036673 is not accurate 1x : In qbeta(x, p, 1/2) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.0007363 is not accurate 1x : In qbeta(x, p, 1/2) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.014403 is not accurate 1x : In qbeta(x, p, 1/2) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.36491 is not accurate 2x : In qbeta(x, p, 1/2) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.4 is not accurate 1x : In qbeta(x, p, 1/2) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.0073506 is not accurate 1x : In qbeta(x, p, 1/2) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.0053028 is not accurate 1x : In qbeta(x, p, 1/2) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.0051553 is not accurate 1x : In qbeta(x, p, 1/2) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.46491 is not accurate 2x : In qbeta(x, p, 1/2) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.5 is not accurate > showProc.time() Time (user system elapsed): 0.08 0 0.07 > > ### q = p !! > for(x in pr.val) { + cat("-----------\nx=", format01prec(x, digits = 6),":\n~~~~~~~~~~\n") + for(p in p.seq[p.seq >= 1/2]) { #-- otherwise, qt(.) is not defined + y <- 1/(1+ qt(1-x, 2*p)^2/(2*p)) + qb0 <- 1/2 + ifelse(x < 1/2, -1, 1)* sqrt(1-y)/2 #-- should be = qbeta(..) + t0 <- proc.time()[1] + + qb <- qbeta(x, p,p) + C1 <- (t1 <- proc.time()[1]) - t0 + pb <- pbeta(qb, p,p) + C2 <- (t2 <- proc.time()[1]) - t1 + ## NO Debug: + fcat(p,paste("|",t2-t0,"|",sep = ''),format01prec(qb,digits = 10), + paste("|",formatC(qb0-qb,dig = 8,w = 12)), qb/qb0 - 1, pb/x - 1, f.wid = 12) + ## Debug: + ##fcat(format01prec(qb,digits=10,flag='-'),paste("|",formatC(qb0-qb,dig=8,w=12)), + ## qb/qb0 - 1, pb/x - 1, f.wid=12) + } + } ----------- x= 1e-05 : ~~~~~~~~~~ 1 |0| 1e-05 | -4.551108e-17 4.551e-12 -3.331e-16 10 |0| 0.1108656971 | -5.6038507e-14 5.054e-13 -4.108e-15 100 |0| 0.3524028603 | -3.3584246e-14 9.526e-14 -3.553e-15 1e+10 |0| 0.4999849213 | -4.306e-13 8.613e-13 -3.141e-11 1e+11 |0| 0.4999952317 | 2.7752245e-12 -5.55e-12 -7.884e-11 1e+12 |0| 0.4999984921 | -6.4852568e-12 1.297e-11 -1.241e-10 1e+13 |0| 0.4999995232 | -7.8716478e-12 1.574e-11 2.711e-09 1e+14 |0| 0.4999998492 | -7.6098794e-11 1.522e-10 -1.011e-09 1e+15 |0| 0.4999999523 | -2.4064528e-11 4.813e-11 -1.137e-08 1e+16 |0| 0.4999999849 | 1.775049e-10 -3.55e-10 -1.178e-07 1e+17 |0| 1-NaN | NaN NaN NaN 1e+18 |0| 0.4999999985 | 1.5088603e-09 -3.018e-09 -0.01251 ----------- x= 0.0001 : ~~~~~~~~~~ 1 |0| 0.0001 | -1.1018205e-17 1.101e-13 8.882e-16 10 |0| 0.143851514 | -1.7486013e-15 1.221e-14 0 100 |0| 0.3705985398 | -9.9920072e-16 2.665e-15 -9.992e-16 1e+10 |0| 0.4999868513 | 3.1796787e-13 -6.359e-13 4.908e-11 1e+11 |0| 0.499995842 | 2.1037061e-12 -4.207e-12 2.365e-10 1e+12 |0| 0.4999986851 | -3.5565439e-12 7.113e-12 4.647e-10 1e+13 |0| 0.4999995842 | 3.2252812e-11 -6.451e-11 1.767e-09 1e+14 |0| 0.4999998685 | 9.4668828e-11 -1.893e-10 -1.972e-09 1e+15 |0| 0.4999999584 | 9.6791242e-11 -1.936e-10 -7.221e-09 1e+16 |0| 0.4999999869 | 2.4392455e-10 -4.878e-10 2.027e-07 1e+17 |0| 1-NaN | NaN NaN NaN 1e+18 |0| 0.4999999987 | 1.3158474e-09 -2.632e-09 -0.01088 ----------- x= 0.001 : ~~~~~~~~~~ 1 |0| 0.001 | 8.6736174e-19 -8.882e-16 2.22e-16 10 |0| 0.1890371159 | 5.5511151e-17 -3.331e-16 1.998e-15 100 |0| 0.3919037825 | -5.5511151e-17 2.22e-16 1.998e-15 1e+10 |0| 0.4999890744 | 6.9078077e-13 -1.382e-12 -1.479e-11 1e+11 |0| 0.499996545 | 3.4321435e-12 -6.864e-12 1.013e-10 1e+12 |0| 0.4999989074 | -9.0762398e-12 1.815e-11 6.884e-11 1e+13 |0| 0.4999996545 | 2.9265035e-11 -5.853e-11 -1.581e-09 1e+14 |0| 0.4999998907 | 9.2544306e-12 -1.851e-11 1.824e-09 1e+15 |0| 0.4999999655 | -3.9647319e-10 7.929e-10 4.302e-09 1e+16 |0| 0.4999999891 | 3.8890896e-10 -7.778e-10 -4.128e-09 1e+17 |0| 0.4999999965 | 3.4549842e-09 -6.91e-09 1.781e-06 1e+18 |0| 0.4999999989 | 1.0934074e-09 -2.187e-09 -0.00802 ----------- x= 0.01 : ~~~~~~~~~~ 1 |0| 0.01 | 8.6736174e-18 -8.882e-16 4.441e-16 10 |0| 0.2539530795 | -5.5511151e-17 2.22e-16 2.22e-16 100 |0| 0.4182038851 | 1.110223e-16 -2.22e-16 -8.882e-16 1e+10 |0| 0.4999917751 | 9.6855857e-13 -1.937e-12 7.777e-12 1e+11 |0| 0.4999973991 | 3.0642155e-13 -6.128e-13 1.494e-11 1e+12 |0| 0.4999991775 | -1.6775858e-11 3.355e-11 -1.713e-10 1e+13 |0| 0.4999997399 | -3.7316539e-11 7.463e-11 -7.878e-10 1e+14 |0| 0.4999999178 | -4.5534687e-11 9.107e-11 3.329e-09 1e+15 |0| 0.499999974 | 1.9979163e-10 -3.996e-10 1.271e-08 1e+16 |0| 0.4999999918 | 7.7430118e-10 -1.549e-09 -1.081e-08 1e+17 |0| 1-NaN | NaN NaN NaN 1e+18 |0| 0.4999999992 | 8.229486e-10 -1.646e-09 -0.003466 ----------- x= 0.1 : ~~~~~~~~~~ 1 |0| 0.1 | -2.7755576e-17 2.22e-16 2.22e-16 10 |0| 0.3579298802 | 5.5511151e-17 -1.11e-16 -7.772e-16 100 |0| 0.4547268852 | -1.110223e-16 2.22e-16 1.11e-15 1e+10 |0| 0.499995469 | -9.7882813e-13 1.958e-12 -8.605e-12 1e+11 |0| 0.4999985672 | -3.0953018e-13 6.191e-13 -7.554e-12 1e+12 |0| 0.4999995469 | 1.8279711e-11 -3.656e-11 -1.327e-10 1e+13 |0| 0.4999998567 | -3.2958469e-11 6.592e-11 9.523e-10 1e+14 |0| 0.4999999547 | -1.0422385e-11 2.084e-11 2.17e-09 1e+15 |0| 0.4999999857 | -5.7297911e-10 1.146e-09 3.712e-09 1e+16 |0| 0.4999999955 | 4.5309689e-09 -9.062e-09 5.569e-08 1e+17 |0| 1-NaN | NaN NaN NaN 1e+18 |0| 0.4987512411 | 0.0012487589 -0.002498 -1 ----------- x= 0.2 : ~~~~~~~~~~ 1 |0| 0.2 | 0 0 0 10 |0| 0.4055828363 | -5.5511151e-17 2.22e-16 -5.551e-16 100 |0| 0.4702334145 | 4.4408921e-16 -9.992e-16 -2.22e-15 1e+10 |0| 0.4999970244 | -9.0277785e-13 1.806e-12 9.034e-12 1e+11 |0| 0.499999059 | 2.6642022e-12 -5.328e-12 -2.721e-12 1e+12 |0| 0.4999997024 | 8.4249274e-13 -1.685e-12 5.764e-11 1e+13 |0| 0.4999999059 | -1.4710383e-10 2.942e-10 3.216e-10 1e+14 |0| 0.4999999702 | -4.6518289e-11 9.304e-11 3.959e-10 1e+15 |0| 0.4999999906 | -1.1271006e-09 2.254e-09 -6.525e-09 1e+16 |0| 0.499999997 | 2.9755804e-09 -5.951e-09 2.87e-09 1e+17 |0| 1-NaN | NaN NaN NaN 1e+18 |0| 0.4999999997 | 2.9666092e-10 -5.933e-10 0.003556 ----------- x= 0.4 : ~~~~~~~~~~ 1 |0| 0.4 | -5.5511151e-17 2.22e-16 0 10 |0| 0.471342474 | -5.5511151e-17 2.22e-16 4.441e-16 100 |0| 0.4910323454 | -5.5511151e-16 1.11e-15 4.441e-16 1e+10 |0| 0.4999991043 | 3.7760906e-12 -7.552e-12 -6.288e-12 1e+11 |0| 0.4999997167 | 3.059264e-11 -6.119e-11 -3.867e-11 1e+12 |0| 0.4999999104 | -1.4514362e-10 2.903e-10 4.298e-11 1e+13 |0| 0.4999999717 | 4.4754728e-10 -8.951e-10 3.301e-10 1e+14 |0| 0.499999991 | 1.5065921e-09 -3.013e-09 -9.044e-11 1e+15 |0| 0.4999999972 | 2.8325067e-09 -5.665e-09 -2.37e-09 1e+16 |0| 0.4999999991 | 8.9571706e-10 -1.791e-09 4.729e-08 1e+17 |0| 1-NaN | NaN NaN NaN 1e+18 |0| 0.4999999999 | 8.884482e-11 -1.777e-10 0.001986 ----------- x= 0.5 : ~~~~~~~~~~ 1 |0| 0.5 | 0 0 0 10 |0| 0.5 | 5.5511151e-17 -1.11e-16 -7.772e-16 100 |0| 0.5 | 0 0 1.554e-15 1e+10 |0| 0.5 | 0 0 0 1e+11 |0| 0.5 | 0 0 0 1e+12 |0| 0.5 | 0 0 2.22e-16 1e+13 |0| 0.5 | 0 0 0 1e+14 |0| 0.5 | 0 0 2.22e-16 1e+15 |0| 0.5 | 1.110223e-16 -2.22e-16 -8.921e-09 1e+16 |0| 0.5 | -4.4408921e-16 8.882e-16 9.027e-08 1e+17 |0| 0.5 | -1.110223e-16 2.22e-16 5.709e-08 1e+18 |0| 0.5 | 1.2212453e-15 -2.442e-15 -2.744e-06 ----------- x= 0.6 : ~~~~~~~~~~ 1 |0| 0.6 | 1.110223e-16 -2.22e-16 0 10 |0| 0.528657526 | 1.110223e-16 -2.22e-16 -2.22e-16 100 |0| 0.5089676546 | 6.6613381e-16 -1.332e-15 -1.665e-15 1e+10 |0| 0.5000008957 | -3.7760906e-12 7.552e-12 4.192e-12 1e+11 |0| 0.5000002833 | -3.0592529e-11 6.119e-11 -1.816e-11 1e+12 |0| 0.5000000896 | 1.4514356e-10 -2.903e-10 8.25e-11 1e+13 |0| 0.5000000283 | -4.4754722e-10 8.951e-10 -2.201e-10 1e+14 |0| 0.500000009 | -1.5065922e-09 3.013e-09 1.483e-09 1e+15 |0| 0.5000000028 | -2.8325067e-09 5.665e-09 -2.02e-09 1e+16 |0| 0.5000000009 | -8.9571706e-10 1.791e-09 -3.153e-08 1e+17 |0| 1-NaN | NaN NaN NaN 1e+18 |0| 0.5000000001 | -8.884482e-11 1.777e-10 -0.001324 ----------- x= 0.8 : ~~~~~~~~~~ 1 |0| 0.8 | 0 0 0 10 |0| 0.5944171637 | 0 0 0 100 |0| 0.5297665855 | -4.4408921e-16 8.882e-16 4.441e-16 1e+10 |0| 0.5000029756 | 9.0283336e-13 -1.806e-12 -1.17e-11 1e+11 |0| 0.500000941 | -2.6642022e-12 5.328e-12 6.803e-13 1e+12 |0| 0.5000002976 | -8.4243723e-13 1.685e-12 -7.482e-11 1e+13 |0| 0.5000000941 | 1.4710388e-10 -2.942e-10 -8.041e-11 1e+14 |0| 0.5000000298 | 4.6518345e-11 -9.304e-11 -9.898e-11 1e+15 |0| 0.5000000094 | 1.1271006e-09 -2.254e-09 3.588e-09 1e+16 |0| 0.500000003 | -2.9755804e-09 5.951e-09 -7.176e-10 1e+17 |0| 1-NaN | NaN NaN NaN 1e+18 |0| 0.5000000003 | -2.9666092e-10 5.933e-10 -0.0008889 ----------- x= 0.9 : ~~~~~~~~~~ 1 |0| 0.9 | 0 0 0 10 |0| 0.6420701198 | 0 0 2.22e-16 100 |0| 0.5452731148 | 0 0 4.441e-16 1e+10 |0| 0.500004531 | 9.7888364e-13 -1.958e-12 9.561e-13 1e+11 |0| 0.5000014328 | 3.0953018e-13 -6.191e-13 8.393e-13 1e+12 |0| 0.5000004531 | -1.8279711e-11 3.656e-11 1.475e-11 1e+13 |0| 0.5000001433 | 3.2958414e-11 -6.592e-11 6.451e-11 1e+14 |0| 0.5000000453 | 1.042233e-11 -2.084e-11 -2.411e-10 1e+15 |0| 0.5000000143 | 5.7297911e-10 -1.146e-09 -4.124e-10 1e+16 |0| 0.5000000045 | -4.5309689e-09 9.062e-09 -6.187e-09 1e+17 |0| 1-NaN | NaN NaN NaN 1e+18 |0| 0.5012487589 | -0.0012487589 0.002498 0.1111 ----------- x= 0.99 : ~~~~~~~~~~ 1 |0| 0.99 | 0 0 0 10 |0| 0.7460469205 | 0 0 0 100 |0| 0.5817961149 | -1.110223e-16 2.22e-16 0 1e+10 |0| 0.5000082249 | -9.6855857e-13 1.937e-12 -7.849e-14 1e+11 |0| 0.5000026009 | -3.0642155e-13 6.128e-13 -1.509e-13 1e+12 |0| 0.5000008225 | 1.6775914e-11 -3.355e-11 1.731e-12 1e+13 |0| 0.5000002601 | 3.7316483e-11 -7.463e-11 7.958e-12 1e+14 |0| 0.5000000822 | 4.5534687e-11 -9.107e-11 -9.311e-11 1e+15 |0| 0.500000026 | -1.9979163e-10 3.996e-10 -1.284e-10 1e+16 |0| 0.5000000082 | -7.7430118e-10 1.549e-09 1.091e-10 1e+17 |0| 1-NaN | NaN NaN NaN 1e+18 |0| 0.5000000008 | -8.229486e-10 1.646e-09 3.501e-05 ----------- x= 0.999 : ~~~~~~~~~~ 1 |0| 0.999 | 0 0 0 10 |0| 0.8109628841 | -1.110223e-16 2.22e-16 0 100 |0| 0.6080962175 | 0 0 0 1e+10 |0| 0.5000109256 | -6.9078077e-13 1.382e-12 1.488e-14 1e+11 |0| 0.500003455 | -3.4321435e-12 6.864e-12 1.286e-13 1e+12 |0| 0.5000010926 | 9.0762953e-12 -1.815e-11 -6.894e-14 1e+13 |0| 0.5000003455 | -2.9265035e-11 5.853e-11 1.583e-12 1e+14 |0| 0.5000001093 | -9.254375e-12 1.851e-11 -1.825e-12 1e+15 |0| 0.5000000345 | 3.9647319e-10 -7.929e-10 -2.315e-11 1e+16 |0| 0.5000000109 | -3.8890891e-10 7.778e-10 4.133e-12 1e+17 |0| 0.5000000035 | -3.4549842e-09 6.91e-09 -1.783e-09 1e+18 |0| 0.5000000011 | -1.0934074e-09 2.187e-09 8.028e-06 ----------- x= 0.9999 : ~~~~~~~~~~ 1 |0| 0.9999 | 0 0 0 10 |0| 0.856148486 | -1.110223e-16 2.22e-16 0 100 |0| 0.6294014602 | 1.110223e-16 -2.22e-16 0 1e+10 |0| 0.5000131487 | -3.1796787e-13 6.359e-13 -4.885e-15 1e+11 |0| 0.500004158 | -2.1036506e-12 4.207e-12 -5.063e-14 1e+12 |0| 0.5000013149 | 3.5565995e-12 -7.113e-12 -4.641e-14 1e+13 |0| 0.5000004158 | -3.2252756e-11 6.451e-11 -1.767e-13 1e+14 |0| 0.5000001315 | -9.4668828e-11 1.893e-10 -6.776e-13 1e+15 |0| 0.5000000416 | -9.6791353e-11 1.936e-10 2.935e-12 1e+16 |0| 0.5000000131 | -2.4392455e-10 4.878e-10 -2.027e-11 1e+17 |0| 1-NaN | NaN NaN NaN 1e+18 |0| 0.5000000013 | -1.3158474e-09 2.632e-09 1.088e-06 ----------- x= 0.99999 : ~~~~~~~~~~ 1 |0| 0.99999 | 0 0 0 10 |0| 0.8891343029 | 0 0 0 100 |0| 0.6475971397 | 0 0 0 1e+10 |0| 0.5000150787 | 4.3054449e-13 -8.611e-13 4.441e-16 1e+11 |0| 0.5000047683 | -2.7752245e-12 5.55e-12 3.775e-15 1e+12 |0| 0.5000015079 | 6.4852568e-12 -1.297e-11 1.332e-15 1e+13 |0| 0.5000004768 | 7.8717033e-12 -1.574e-11 -2.72e-14 1e+14 |0| 0.5000001508 | 7.6098794e-11 -1.522e-10 1.021e-14 1e+15 |0| 0.5000000477 | 2.4064528e-11 -4.813e-11 1.137e-13 1e+16 |0| 0.5000000151 | -1.775049e-10 3.55e-10 1.178e-12 1e+17 |0| 1-NaN | NaN NaN NaN 1e+18 |0| 0.5000000015 | -1.5088603e-09 3.018e-09 1.251e-07 There were 18 warnings (use warnings() to see them) > > summary(warnings()) Summary of (a total of 18) warning messages: 12x : In qbeta(x, p, p) : NaNs produced 2x : In qbeta(x, p, p) : qbeta(a, *) =: x0 with |pbeta(x0,*, log_) - alpha| = 6.2376e+12 is not accurate 2x : In qbeta(x, p, p) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.00071114 is not accurate 2x : In qbeta(x, p, p) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.00079455 is not accurate > showProc.time() Time (user system elapsed): 0.07 0.02 0.08 > > > ###---- 2008-11-12 : another example, not sure if covered above: > ### Letting a --> 0 --- Distribution only defined for a > 0 --- > ## > curve(qbeta(.95, x, 20), 1e-300, 0.1, n=1001)# a -> 0 : qbeta() -> 0 > > ## now, "zoom in logarithmically": > curve(qbeta(.95, x, 20), 1e-7, 0.1, log="xy", n=1001, yaxt="n") > if(require("sfsmisc")) { eaxis(2, at= c(10^-c(314,309), axTicks(2, nintLog=15))) + } else { axis(2) + warning("could not use the cool eaxis() function from CRAN package 'sfsmisc' as that is not yet installed.") + } Loading required package: sfsmisc Attaching package: 'sfsmisc' The following object is masked _by_ '.GlobalEnv': list_ > > ## -> clearly should go to 0; > ## now does ... not continually, as with pbeta() warning about bgrat() underflow > ## now without any warnings and continually: not 0 but staying at > unique(qbeta(.95, 10^-(100:5), 20)) [1] 5.562685e-309 > ## 5.5626846462680034577e-309 > > ## limit : > qbeta(.95, 0, 20)# gave '1' instead of 0 --> now 0 {as for dbeta(*, a=0, *) [1] 0 > > ## Variation on the above (a --> 0) : > ## --------- > ## this seems ok; is "the" behavior at large" > x. <- seq(0,1,len=1001); plot(x., pbeta(x., 1e-20,20), type="l") > ## but if you log-zoom in , > ## the precision problem is already in pbeta(): ----- fixed on 2009-10-16 (in Skafidia, Greece) > ##-- TOMS 708: 'a0 small .. -> apser() > curve(pbeta(exp(x), 1e-16, 10, log.p=TRUE), -10, 0, n=1000) > curve(pbeta(exp(x), 1e-16, 10, log.p=TRUE), -900, 0, n=1000) > ## _jumps_ to 0 - because exp() *must* [no bug in pbeta !] > ## but this is somewhat similar, *not* using apser(), but L20, then > ## a) (x < 0.29): bup() + L131 bgrat() & w = 1-w1 ~= 1 -- now fixed > ## b) x >=0.29: bpser() > xpb <- curve(pbeta(exp(x), 2e-15, 10, log.p=TRUE), -21, 0, n=1000) > xpb <- within(as.data.frame(xpb), {exp.x <- exp(x); bp <- exp(x) >= 0.29}) > ## unfinished > bps <- bpser(2e-15, 10, x = with(xpb, exp.x[bp]), log.p=TRUE, verbose=TRUE) bpser(a=2e-15, b=10, x=0.295461, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 3.2194980694981e-15; -> n=4 iterations, |w|=0.240054 <= 0.5=tol:=eps/a ==> a*sum=-3.13989e-15 > -1 bpser(a=2e-15, b=10, x=0.301737, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 3.2615401115401e-15; -> n=4 iterations, |w|=0.261112 <= 0.5=tol:=eps/a ==> a*sum=-3.16983e-15 > -1 bpser(a=2e-15, b=10, x=0.308147, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 3.3035821535822e-15; -> n=4 iterations, |w|=0.284017 <= 0.5=tol:=eps/a ==> a*sum=-3.19881e-15 > -1 bpser(a=2e-15, b=10, x=0.314693, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 3.3456241956242e-15; -> n=4 iterations, |w|=0.308931 <= 0.5=tol:=eps/a ==> a*sum=-3.22669e-15 > -1 bpser(a=2e-15, b=10, x=0.321379, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 3.3876662376662e-15; -> n=4 iterations, |w|=0.336031 <= 0.5=tol:=eps/a ==> a*sum=-3.25335e-15 > -1 bpser(a=2e-15, b=10, x=0.328206, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 3.4297082797083e-15; -> n=4 iterations, |w|=0.365507 <= 0.5=tol:=eps/a ==> a*sum=-3.27863e-15 > -1 bpser(a=2e-15, b=10, x=0.335178, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 3.4717503217503e-15; -> n=4 iterations, |w|=0.39757 <= 0.5=tol:=eps/a ==> a*sum=-3.30237e-15 > -1 bpser(a=2e-15, b=10, x=0.342299, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 3.5137923637924e-15; -> n=4 iterations, |w|=0.432445 <= 0.5=tol:=eps/a ==> a*sum=-3.32439e-15 > -1 bpser(a=2e-15, b=10, x=0.34957, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 3.5558344058344e-15; -> n=4 iterations, |w|=0.470379 <= 0.5=tol:=eps/a ==> a*sum=-3.34449e-15 > -1 bpser(a=2e-15, b=10, x=0.356996, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 3.5978764478764e-15; -> n=5 iterations, |w|=0.146123 <= 0.5=tol:=eps/a ==> a*sum=-3.65471e-15 > -1 bpser(a=2e-15, b=10, x=0.36458, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 3.6399184899185e-15; -> n=5 iterations, |w|=0.162317 <= 0.5=tol:=eps/a ==> a*sum=-3.70269e-15 > -1 bpser(a=2e-15, b=10, x=0.372325, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 3.6819605319605e-15; -> n=5 iterations, |w|=0.180307 <= 0.5=tol:=eps/a ==> a*sum=-3.75163e-15 > -1 bpser(a=2e-15, b=10, x=0.380235, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 3.7240025740026e-15; -> n=5 iterations, |w|=0.20029 <= 0.5=tol:=eps/a ==> a*sum=-3.80163e-15 > -1 bpser(a=2e-15, b=10, x=0.388312, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 3.7660446160446e-15; -> n=5 iterations, |w|=0.222487 <= 0.5=tol:=eps/a ==> a*sum=-3.85281e-15 > -1 bpser(a=2e-15, b=10, x=0.396561, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 3.8080866580867e-15; -> n=5 iterations, |w|=0.247145 <= 0.5=tol:=eps/a ==> a*sum=-3.90531e-15 > -1 bpser(a=2e-15, b=10, x=0.404986, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 3.8501287001287e-15; -> n=5 iterations, |w|=0.274535 <= 0.5=tol:=eps/a ==> a*sum=-3.9593e-15 > -1 bpser(a=2e-15, b=10, x=0.413589, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 3.8921707421707e-15; -> n=5 iterations, |w|=0.304961 <= 0.5=tol:=eps/a ==> a*sum=-4.01495e-15 > -1 bpser(a=2e-15, b=10, x=0.422375, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 3.9342127842128e-15; -> n=5 iterations, |w|=0.338759 <= 0.5=tol:=eps/a ==> a*sum=-4.07247e-15 > -1 bpser(a=2e-15, b=10, x=0.431348, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 3.9762548262548e-15; -> n=5 iterations, |w|=0.376303 <= 0.5=tol:=eps/a ==> a*sum=-4.13209e-15 > -1 bpser(a=2e-15, b=10, x=0.440511, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 4.0182968682969e-15; -> n=5 iterations, |w|=0.418008 <= 0.5=tol:=eps/a ==> a*sum=-4.19407e-15 > -1 bpser(a=2e-15, b=10, x=0.449869, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 4.0603389103389e-15; -> n=5 iterations, |w|=0.464335 <= 0.5=tol:=eps/a ==> a*sum=-4.25871e-15 > -1 bpser(a=2e-15, b=10, x=0.459426, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 4.102380952381e-15; -> n=6 iterations, |w|=0.13165 <= 0.5=tol:=eps/a ==> a*sum=-4.06305e-15 > -1 bpser(a=2e-15, b=10, x=0.469186, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 4.144422994423e-15; -> n=6 iterations, |w|=0.149347 <= 0.5=tol:=eps/a ==> a*sum=-4.09867e-15 > -1 bpser(a=2e-15, b=10, x=0.479153, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 4.186465036465e-15; -> n=6 iterations, |w|=0.169423 <= 0.5=tol:=eps/a ==> a*sum=-4.13334e-15 > -1 bpser(a=2e-15, b=10, x=0.489332, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 4.2285070785071e-15; -> n=6 iterations, |w|=0.192198 <= 0.5=tol:=eps/a ==> a*sum=-4.1669e-15 > -1 bpser(a=2e-15, b=10, x=0.499727, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 4.2705491205491e-15; -> n=6 iterations, |w|=0.218034 <= 0.5=tol:=eps/a ==> a*sum=-4.19918e-15 > -1 bpser(a=2e-15, b=10, x=0.510343, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 4.3125911625912e-15; -> n=6 iterations, |w|=0.247343 <= 0.5=tol:=eps/a ==> a*sum=-4.22997e-15 > -1 bpser(a=2e-15, b=10, x=0.521184, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 4.3546332046332e-15; -> n=6 iterations, |w|=0.280592 <= 0.5=tol:=eps/a ==> a*sum=-4.25904e-15 > -1 bpser(a=2e-15, b=10, x=0.532256, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 4.3966752466752e-15; -> n=6 iterations, |w|=0.318311 <= 0.5=tol:=eps/a ==> a*sum=-4.28611e-15 > -1 bpser(a=2e-15, b=10, x=0.543563, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 4.4387172887173e-15; -> n=6 iterations, |w|=0.3611 <= 0.5=tol:=eps/a ==> a*sum=-4.31087e-15 > -1 bpser(a=2e-15, b=10, x=0.55511, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 4.4807593307593e-15; -> n=6 iterations, |w|=0.409641 <= 0.5=tol:=eps/a ==> a*sum=-4.33297e-15 > -1 bpser(a=2e-15, b=10, x=0.566903, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 4.5228013728014e-15; -> n=6 iterations, |w|=0.464707 <= 0.5=tol:=eps/a ==> a*sum=-4.35199e-15 > -1 bpser(a=2e-15, b=10, x=0.578946, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 4.5648434148434e-15; -> n=7 iterations, |w|=0.112116 <= 0.5=tol:=eps/a ==> a*sum=-4.59167e-15 > -1 bpser(a=2e-15, b=10, x=0.591245, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 4.6068854568855e-15; -> n=7 iterations, |w|=0.12989 <= 0.5=tol:=eps/a ==> a*sum=-4.63856e-15 > -1 bpser(a=2e-15, b=10, x=0.603805, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 4.6489274989275e-15; -> n=7 iterations, |w|=0.15048 <= 0.5=tol:=eps/a ==> a*sum=-4.68634e-15 > -1 bpser(a=2e-15, b=10, x=0.616632, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 4.6909695409695e-15; -> n=7 iterations, |w|=0.174335 <= 0.5=tol:=eps/a ==> a*sum=-4.73516e-15 > -1 bpser(a=2e-15, b=10, x=0.629731, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 4.7330115830116e-15; -> n=7 iterations, |w|=0.201972 <= 0.5=tol:=eps/a ==> a*sum=-4.78521e-15 > -1 bpser(a=2e-15, b=10, x=0.643109, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 4.7750536250536e-15; -> n=7 iterations, |w|=0.233989 <= 0.5=tol:=eps/a ==> a*sum=-4.83672e-15 > -1 bpser(a=2e-15, b=10, x=0.656771, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 4.8170956670957e-15; -> n=7 iterations, |w|=0.271082 <= 0.5=tol:=eps/a ==> a*sum=-4.88994e-15 > -1 bpser(a=2e-15, b=10, x=0.670723, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 4.8591377091377e-15; -> n=7 iterations, |w|=0.314055 <= 0.5=tol:=eps/a ==> a*sum=-4.94519e-15 > -1 bpser(a=2e-15, b=10, x=0.684971, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 4.9011797511798e-15; -> n=7 iterations, |w|=0.36384 <= 0.5=tol:=eps/a ==> a*sum=-5.00284e-15 > -1 bpser(a=2e-15, b=10, x=0.699522, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 4.9432217932218e-15; -> n=7 iterations, |w|=0.421518 <= 0.5=tol:=eps/a ==> a*sum=-5.06331e-15 > -1 bpser(a=2e-15, b=10, x=0.714383, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 4.9852638352638e-15; -> n=7 iterations, |w|=0.488339 <= 0.5=tol:=eps/a ==> a*sum=-5.12712e-15 > -1 bpser(a=2e-15, b=10, x=0.729559, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 5.0273058773059e-15; -> n=8 iterations, |w|=0.090289 <= 0.5=tol:=eps/a ==> a*sum=-5.01429e-15 > -1 bpser(a=2e-15, b=10, x=0.745057, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 5.0693479193479e-15; -> n=8 iterations, |w|=0.106824 <= 0.5=tol:=eps/a ==> a*sum=-5.05363e-15 > -1 bpser(a=2e-15, b=10, x=0.760885, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 5.11138996139e-15; -> n=8 iterations, |w|=0.126387 <= 0.5=tol:=eps/a ==> a*sum=-5.09239e-15 > -1 bpser(a=2e-15, b=10, x=0.777049, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 5.153432003432e-15; -> n=8 iterations, |w|=0.149533 <= 0.5=tol:=eps/a ==> a*sum=-5.13048e-15 > -1 bpser(a=2e-15, b=10, x=0.793556, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 5.195474045474e-15; -> n=8 iterations, |w|=0.176918 <= 0.5=tol:=eps/a ==> a*sum=-5.16774e-15 > -1 bpser(a=2e-15, b=10, x=0.810414, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 5.2375160875161e-15; -> n=8 iterations, |w|=0.209318 <= 0.5=tol:=eps/a ==> a*sum=-5.20401e-15 > -1 bpser(a=2e-15, b=10, x=0.82763, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 5.2795581295581e-15; -> n=8 iterations, |w|=0.247652 <= 0.5=tol:=eps/a ==> a*sum=-5.23907e-15 > -1 bpser(a=2e-15, b=10, x=0.845212, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 5.3216001716002e-15; -> n=8 iterations, |w|=0.293006 <= 0.5=tol:=eps/a ==> a*sum=-5.27268e-15 > -1 bpser(a=2e-15, b=10, x=0.863167, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 5.3636422136422e-15; -> n=8 iterations, |w|=0.346666 <= 0.5=tol:=eps/a ==> a*sum=-5.30453e-15 > -1 bpser(a=2e-15, b=10, x=0.881504, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 5.4056842556843e-15; -> n=8 iterations, |w|=0.410153 <= 0.5=tol:=eps/a ==> a*sum=-5.33427e-15 > -1 bpser(a=2e-15, b=10, x=0.90023, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 5.4477262977263e-15; -> n=8 iterations, |w|=0.485266 <= 0.5=tol:=eps/a ==> a*sum=-5.36143e-15 > -1 bpser(a=2e-15, b=10, x=0.919354, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 5.4897683397683e-15; -> n=9 iterations, |w|=0.0521317 <= 0.5=tol:=eps/a ==> a*sum=-5.48977e-15 > -1 bpser(a=2e-15, b=10, x=0.938884, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 5.5318103818104e-15; -> n=9 iterations, |w|=0.0629892 <= 0.5=tol:=eps/a ==> a*sum=-5.53181e-15 > -1 bpser(a=2e-15, b=10, x=0.958829, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 5.5738524238524e-15; -> n=9 iterations, |w|=0.0761079 <= 0.5=tol:=eps/a ==> a*sum=-5.57385e-15 > -1 bpser(a=2e-15, b=10, x=0.979198, log=1, eps=1e-15): log(x^a/(a*B(a,b))) = 5.6158944658945e-15; -> n=9 iterations, |w|=0.0919589 <= 0.5=tol:=eps/a ==> a*sum=-5.61589e-15 > -1 > > > u <- seq(-16, 1.25, length=501) > op <- mult.fig(mfrow=c(3,5), marP = c(0,-1.5,-1,-.8), + main = expression(pbeta(e^u, alpha, beta, ~~ log=='TRUE')))$old.par > for(b in c(.5, 2, 5)) + for(a in c(.1, .2, .5, 1, 2)/100) { + plot(u, (qp <- pbeta(exp(u), a,b, log.p=TRUE)), ylab = NA, + type="l", main = substitute(list(alpha == a, beta == b), list(a=a,b=b))) + } > par(op) > > > summary(warnings()) Summary of (a total of 18) warning messages: 12x : In qbeta(x, p, p) : NaNs produced 2x : In qbeta(x, p, p) : qbeta(a, *) =: x0 with |pbeta(x0,*, log_) - alpha| = 6.2376e+12 is not accurate 2x : In qbeta(x, p, p) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.00071114 is not accurate 2x : In qbeta(x, p, p) : qbeta(a, *) =: x0 with |pbeta(x0,*) - alpha| = 0.00079455 is not accurate > showProc.time() Time (user system elapsed): 0.06 0 0.06 > > ###---- Another (or the same?) set of problems: > ## From: Martin Maechler > ## Sender: r-devel-bounces@r-project.org > ## To: leydold@statistik.wu-wien.ac.at > ## Cc: r-devel@r-project.org > ## Subject: Re: [Rd] Buglet in qbeta? > ## Date: Wed, 7 Oct 2009 16:04:47 +0200 > > ## >>>>> "JL" == Josef Leydold > ## >>>>> on Wed, 7 Oct 2009 14:48:26 +0200 writes: > > ## JL> Hi, > ## JL> I sometimes play around with extreme parameters for distributions and > ## JL> found that qbeta is not always monotone as the following example shows. > ## JL> I don't know whether this is serious enough to submit a bug report (as > ## JL> this example is near to the limitations of floating point arithmetic). > > ## well, it's not near enough the limits to *not* warrant a bug > ## report! > > ## I "know" (I have saved a collection of ..) about several > ## accuracy problems of pbeta() / qbeta(), most cases indeed > ## "borderline" (at the limits of FP arithmetic) but this one may > ## well be a new one. > > ## A bit more succint than the dump of your numbers is a simple > ## graphic > > p <- (1:100)/100; qp <- qbeta(p, 0.01, 5) > plot(p,qp, type="o", log = "y") # now fine > > ## or, already suggesting a fix to the bug: > > plot(p,qp, type="o", log = "xy")## now fine > ## which is definitely a bug... though not a terrible one. > ## ... > > ## more experiments -- adapt to log-scale: > if(!exists("lseq", mode="function"))# package "sfsmisc" + lseq <- function (from, to, length) exp(seq(log(from), log(to), length.out = length)) > p <- lseq(.005,1, len = 500) ; qp <- qbeta(p, 0.01,5) > ## (no warnings) > plot(p,qp, ylab = "qbeta(p, .01, 5)", type="l", log = "xy") > > ## p --> 0 even closer -- next problem {but that's easier forgivable ..} > p <- lseq(.0005,1, len = 500) > a <- .01 ; b <- 5 > ## nomore (gave 36 warnings "qbeta: ..precision not reached ..") > plot(p, (qp <- qbeta(p, a,b)), ylab = expression(qbeta(p, alpha, beta)), + type="l", log = "xy", main = substitute(list(alpha == a, beta == b), list(a=a,b=b))) > > ## Seems almost independent of beta =: b > op <- mult.fig(mfrow = c(3,5), marP = c(0,-.8,-1,-.8), + main = quote('qbeta(p, , .) vs p for '~ {p %->% 0}))$old.par > for(b in c(.5, 2, 5)) + for(a in c(.1, .2, .5, 1, 2)/100) { # last one, a = 0.02 does not underflow + plot(p, (qp <- qbeta(p, a,b)), ylab = expression(qbeta(p, alpha, beta)), col=2, + type="l", log = "xy", main = substitute(list(alpha == a, beta == b), list(a=a,b=b))) + } Warning messages: 1: In xy.coords(x, y, xlabel, ylabel, log) : 407 y values <= 0 omitted from logarithmic plot 2: In xy.coords(x, y, xlabel, ylabel, log) : 361 y values <= 0 omitted from logarithmic plot 3: In xy.coords(x, y, xlabel, ylabel, log) : 221 y values <= 0 omitted from logarithmic plot 4: In xy.coords(x, y, xlabel, ylabel, log) : 407 y values <= 0 omitted from logarithmic plot 5: In xy.coords(x, y, xlabel, ylabel, log) : 361 y values <= 0 omitted from logarithmic plot 6: In xy.coords(x, y, xlabel, ylabel, log) : 221 y values <= 0 omitted from logarithmic plot 7: In xy.coords(x, y, xlabel, ylabel, log) : 408 y values <= 0 omitted from logarithmic plot 8: In xy.coords(x, y, xlabel, ylabel, log) : 361 y values <= 0 omitted from logarithmic plot 9: In xy.coords(x, y, xlabel, ylabel, log) : 222 y values <= 0 omitted from logarithmic plot > par(op) > > summary(warnings()) Summary of (a total of 9) warning messages: 2x : In xy.coords(x, y, xlabel, ylabel, log) : 407 y values <= 0 omitted from logarithmic plot 3x : In xy.coords(x, y, xlabel, ylabel, log) : 361 y values <= 0 omitted from logarithmic plot 2x : In xy.coords(x, y, xlabel, ylabel, log) : 221 y values <= 0 omitted from logarithmic plot 1x : In xy.coords(x, y, xlabel, ylabel, log) : 408 y values <= 0 omitted from logarithmic plot 1x : In xy.coords(x, y, xlabel, ylabel, log) : 222 y values <= 0 omitted from logarithmic plot > showProc.time() Time (user system elapsed): 0.06 0 0.07 > > ## x <- qbeta((0:100)/100,0.01,5) > ## x > ## .... > ## order(x) > ## 1 2 3 .... 50 51 55 68 72 56 59 ... > ## pbeta(x,0.01,5) > ## ... > ## >> version > ## JL> _ > ## JL> platform x86_64-unknown-linux-gnu > ## JL> arch x86_64 > ## JL> os linux-gnu > ## JL> system x86_64, linux-gnu > ## JL> status Under development (unstable) > ## JL> major 2 > ## JL> minor 11.0 > ## JL> year 2009 > ## JL> month 10 > ## JL> day 07 > ## JL> svn rev 49963 > ## JL> language R > ## JL> version.string R version 2.11.0 Under development (unstable) (2009-10-07 > ## JL> r49963) > > ## JL> p.s. there are similar results for R-2.9.2 in Windows (with > ## JL> different round-off errors). > > > ## --> only vary a == alpha == shape1: > > p.qbeta.p0 <- function(p1 = .005, p2 = 1, beta = 1, mark.p = 0.5, + alphas = c(.1, .15, .2, .3, .4, .5, .8, 1, 2)/100, + n = 500, verbose = getOption("verbose")) + { + ## Purpose: Plot qbeta() investigations for small alpha & p --> 0 + ## ---------------------------------------------------------------------- + ## Arguments: + ## ---------------------------------------------------------------------- + ## Author: Martin Maechler, Date: 7 Oct 2009, 22:46 + stopifnot(exprs = { + require("sfsmisc") # lseq(), mult.fig() + require("DPQ") # qbetaAppr() .. + p1 < p2 + beta > 0 + alphas > 0 ; diff(alphas) > 0 + n == round(n) ; n > 2 + }) + p <- lseq(p1, p2, len = n) + tit <- expression(qbeta(p, alpha, beta) * " for small"~alpha~ + " and " * {p %down% 0}) + cols <- adjustcolor(1:3, 1/2); lwds <- c(2,4,2); ltys <- c(2,1,4) + op <- mult.fig(length(alphas), main = tit, marP = c(-1,-1,-1,-.6))$old.par + on.exit(par(op)) + for(a in alphas) { + matplot(p, cbind(qbetaAppr (p, a,beta, verbose=verbose), + qbeta (p, a,beta), + qbetaAppr.3(p, a,beta)), + col = cols, lty = ltys, lwd = lwds, + ylab = "qbeta", type= "l", log = "xy", + main = substitute(list(alpha == a, beta == b), + list(a=a, b=beta))) + if(match(a, alphas) == length(alphas) %/% 2) + legend("left", c("qbetaAppr", "qbeta", "qbeta.a..3"), + col = cols, lty = ltys, lwd = lwds, inset = .02, bty="n") + if(is.numeric(mark.p) && length(mark.p)) { + abline( v = mark.p, lty = 3, lwd = 2, col = "light blue") + axis(1, at = mark.p, line=-1, lwd = 2, col = "light blue") + } + } + invisible(p) + } > > p.qbeta.p0(.05, mark=NULL, verbose=2) p or q <= 1 : t <= 0: using qbetaAppr.2() t > 0; t' <= 1: qbetaAppr.3() t > 0; t' > 1 : qbetaAppr.4() = 1-2/(t'+1) p or q <= 1 : t <= 0: using qbetaAppr.2() t > 0; t' <= 1: qbetaAppr.3() t > 0; t' > 1 : qbetaAppr.4() = 1-2/(t'+1) p or q <= 1 : t <= 0: using qbetaAppr.2() t > 0; t' <= 1: qbetaAppr.3() t > 0; t' > 1 : qbetaAppr.4() = 1-2/(t'+1) p or q <= 1 : t <= 0: using qbetaAppr.2() t > 0; t' <= 1: qbetaAppr.3() t > 0; t' > 1 : qbetaAppr.4() = 1-2/(t'+1) p or q <= 1 : t <= 0: using qbetaAppr.2() t > 0; t' <= 1: qbetaAppr.3() t > 0; t' > 1 : qbetaAppr.4() = 1-2/(t'+1) p or q <= 1 : t <= 0: using qbetaAppr.2() t > 0; t' <= 1: qbetaAppr.3() t > 0; t' > 1 : qbetaAppr.4() = 1-2/(t'+1) p or q <= 1 : t <= 0: using qbetaAppr.2() t > 0; t' <= 1: qbetaAppr.3() t > 0; t' > 1 : qbetaAppr.4() = 1-2/(t'+1) p or q <= 1 : t <= 0: using qbetaAppr.2() t > 0; t' <= 1: qbetaAppr.3() t > 0; t' > 1 : qbetaAppr.4() = 1-2/(t'+1) p or q <= 1 : t <= 0: using qbetaAppr.2() t > 0; t' <= 1: qbetaAppr.3() t > 0; t' > 1 : qbetaAppr.4() = 1-2/(t'+1) There were 13 warnings (use warnings() to see them) > p.qbeta.p0(.0005) There were 16 warnings (use warnings() to see them) > p.qbeta.p0(.00002) There were 17 warnings (use warnings() to see them) > ## but for beta <~ 0.5, the approximation seems problematic > p.qbeta.p0(beta = 0.8) There were 15 warnings (use warnings() to see them) > p.qbeta.p0(beta = 0.4)# approx. does not go up to 1 There were 27 warnings (use warnings() to see them) > ## but then, the "doc" says to use it for x < 1/2 > p.qbeta.p0(beta = 0.2) There were 27 warnings (use warnings() to see them) > p.qbeta.p0(beta = 0.1)## -> first plot, alpha = .001 is "catastrophe" << still not ok There were 33 warnings (use warnings() to see them) > p.qbeta.p0(.00001, beta = .0001)# here, the approx. is completely hopeless There were 46 warnings (use warnings() to see them) > ## larger alphas : another problem: the *approximation* is bad! > p.qbeta.p0(.001, 0.5, n=2000, alpha= c(.2, .5, 1, 2)) > > showProc.time() Time (user system elapsed): 0.36 0.04 0.4 > > p <- lseq(.0001, 1, len = 2000) > b <- 1 > tit <- expression(qb == qbeta(p, alpha, beta) * " for small"~alpha~" and " * {p %down% 0}) > aset <- c(.1, .15, .2, .3, .4, .5, .8, 1, 2)/100 > op <- mult.fig(length(aset), main = tit, marP = c(-1,-1,-1,-.6))$old.par > c2 <- adjustcolor(2, 1/2) > for(a in aset) { + qpA <- qbetaAppr(p, a,b) + qp <- qbeta (p, a,b) + plot (p, qpA, ylab = "qb", type="l", log = "xy", ylim=pmax(1e-20, range(qp, qpA)), + main = substitute(list(alpha == a, beta == b), list(a=a,b=b))) + lines(p, qp, col = c2, lwd=3) + legend("bottom", c("qbetaAppr()", "qbeta()"), + lty=1, lwd=c(1,3), col=c(palette()[1], c2), bty="n") + } Warning messages: 1: In qbetaAppr(p, a, b) : a[.] > 1/2 (not thought for!) 2: In qbetaAppr(p, a, b) : a[.] > 1/2 (not thought for!) 3: In qbetaAppr(p, a, b) : a[.] > 1/2 (not thought for!) 4: In qbetaAppr(p, a, b) : a[.] > 1/2 (not thought for!) 5: In qbetaAppr(p, a, b) : a[.] > 1/2 (not thought for!) 6: In qbetaAppr(p, a, b) : a[.] > 1/2 (not thought for!) 7: In qbetaAppr(p, a, b) : a[.] > 1/2 (not thought for!) 8: In qbetaAppr(p, a, b) : a[.] > 1/2 (not thought for!) 9: In qbetaAppr(p, a, b) : a[.] > 1/2 (not thought for!) > par(op) > > > ## Ok, let's look at pbeta() "here" > curve(pbeta(x, 0.01, 5, log.p=TRUE), 0, 1, n=1000) # fine (but not close enough to x = 0) > ## log-zoom: > curve(pbeta(x, 0.01, 5, log.p=TRUE), 1e-100, 1, n=10000, log="x") # fine > curve(pbeta(x, 0.01, 5, log.p=TRUE), 1e-250, 1, n=10000, log="x") # fine > ## it does *not* seem to be the problem of pbeta > ## Ok, the "best" qbeta() is the one matching pbeta: > p. <- lseq(1e-3, 1, length=200) > qb <- qbeta (p., 0.02, 0.01) > qb3 <- qbetaAppr.3(p., 0.02, 0.01) > ## qbetaAppr.3() is good for *small* p., but otherwise breaks down *before* > ## ------------ *qbeta()* > cbind(p., pbeta(qb,0.02, 0.01), pbeta(qb3,0.02, 0.01)) p. [1,] 0.001000000 0.001000000 0.3220369 [2,] 0.001035322 0.001035322 0.3222995 [3,] 0.001071891 0.001071891 0.3225640 [4,] 0.001109752 0.001109752 0.3228304 [5,] 0.001148951 0.001148951 0.3230989 [6,] 0.001189534 0.001189534 0.3233695 [7,] 0.001231551 0.001231551 0.3236423 [8,] 0.001275051 0.001275051 0.3239174 [9,] 0.001320088 0.001320088 0.3241949 [10,] 0.001366716 0.001366716 0.3244750 [11,] 0.001414991 0.001414991 0.3247578 [12,] 0.001464971 0.001464971 0.3250435 [13,] 0.001516717 0.001516717 0.3253321 [14,] 0.001570290 0.001570290 0.3256238 [15,] 0.001625756 0.001625756 0.3259188 [16,] 0.001683180 0.001683180 0.3262174 [17,] 0.001742633 0.001742633 0.3265196 [18,] 0.001804186 0.001804186 0.3268258 [19,] 0.001867914 0.001867914 0.3271361 [20,] 0.001933892 0.001933892 0.3274508 [21,] 0.002002200 0.002002200 0.3277701 [22,] 0.002072922 0.002072922 0.3280945 [23,] 0.002146141 0.002146141 0.3284241 [24,] 0.002221947 0.002221947 0.3287594 [25,] 0.002300430 0.002300430 0.3291007 [26,] 0.002381686 0.002381686 0.3294485 [27,] 0.002465811 0.002465811 0.3298033 [28,] 0.002552908 0.002552908 0.3301655 [29,] 0.002643081 0.002643081 0.3305358 [30,] 0.002736440 0.002736440 0.3309148 [31,] 0.002833096 0.002833096 0.3313033 [32,] 0.002933166 0.002933166 0.3317020 [33,] 0.003036771 0.003036771 0.3321120 [34,] 0.003144035 0.003144035 0.3325341 [35,] 0.003255089 0.003255089 0.3329698 [36,] 0.003370064 0.003370064 0.3334202 [37,] 0.003489101 0.003489101 0.3338872 [38,] 0.003612343 0.003612343 0.3343724 [39,] 0.003739937 0.003739937 0.3348780 [40,] 0.003872039 0.003872039 0.3354068 [41,] 0.004008806 0.004008806 0.3359617 [42,] 0.004150405 0.004150405 0.3365466 [43,] 0.004297005 0.004297005 0.3371661 [44,] 0.004448783 0.004448783 0.3378259 [45,] 0.004605922 0.004605922 0.3385333 [46,] 0.004768612 0.004768612 0.3392978 [47,] 0.004937048 0.004937048 0.3401320 [48,] 0.005111433 0.005111433 0.3410530 [49,] 0.005291979 0.005291979 0.3420852 [50,] 0.005478901 0.005478901 0.3432650 [51,] 0.005672426 0.005672426 0.3446498 [52,] 0.005872787 0.005872787 0.3463391 [53,] 0.006080224 0.006080224 0.3485287 [54,] 0.006294989 0.006294989 0.3516974 [55,] 0.006517340 0.006517340 0.3577439 [56,] 0.006747544 0.006747544 1.0000000 [57,] 0.006985880 0.006985880 1.0000000 [58,] 0.007232634 0.007232634 1.0000000 [59,] 0.007488104 0.007488104 1.0000000 [60,] 0.007752597 0.007752597 1.0000000 [61,] 0.008026434 0.008026434 1.0000000 [62,] 0.008309942 0.008309942 1.0000000 [63,] 0.008603464 0.008603464 1.0000000 [64,] 0.008907355 0.008907355 1.0000000 [65,] 0.009221979 0.009221979 1.0000000 [66,] 0.009547716 0.009547716 1.0000000 [67,] 0.009884959 0.009884959 1.0000000 [68,] 0.010234114 0.010234114 1.0000000 [69,] 0.010595602 0.010595602 1.0000000 [70,] 0.010969858 0.010969858 1.0000000 [71,] 0.011357334 0.011357334 1.0000000 [72,] 0.011758496 0.011758496 1.0000000 [73,] 0.012173827 0.012173827 1.0000000 [74,] 0.012603829 0.012603829 1.0000000 [75,] 0.013049020 0.013049020 1.0000000 [76,] 0.013509935 0.013509935 1.0000000 [77,] 0.013987131 0.013987131 1.0000000 [78,] 0.014481182 0.014481182 1.0000000 [79,] 0.014992684 0.014992684 1.0000000 [80,] 0.015522254 0.015522254 1.0000000 [81,] 0.016070528 0.016070528 1.0000000 [82,] 0.016638169 0.016638169 1.0000000 [83,] 0.017225860 0.017225860 1.0000000 [84,] 0.017834309 0.017834309 1.0000000 [85,] 0.018464249 0.018464249 1.0000000 [86,] 0.019116441 0.019116441 1.0000000 [87,] 0.019791669 0.019791669 1.0000000 [88,] 0.020490747 0.020490747 1.0000000 [89,] 0.021214518 0.021214518 1.0000000 [90,] 0.021963854 0.021963854 1.0000000 [91,] 0.022739658 0.022739658 1.0000000 [92,] 0.023542864 0.023542864 1.0000000 [93,] 0.024374442 0.024374442 1.0000000 [94,] 0.025235392 0.025235392 1.0000000 [95,] 0.026126752 0.026126752 1.0000000 [96,] 0.027049597 0.027049597 1.0000000 [97,] 0.028005039 0.028005039 1.0000000 [98,] 0.028994229 0.028994229 1.0000000 [99,] 0.030018358 0.030018358 1.0000000 [100,] 0.031078662 0.031078662 1.0000000 [101,] 0.032176418 0.032176418 1.0000000 [102,] 0.033312948 0.033312948 1.0000000 [103,] 0.034489623 0.034489623 1.0000000 [104,] 0.035707860 0.035707860 1.0000000 [105,] 0.036969127 0.036969127 1.0000000 [106,] 0.038274945 0.038274945 1.0000000 [107,] 0.039626886 0.039626886 1.0000000 [108,] 0.041026581 0.041026581 1.0000000 [109,] 0.042475716 0.042475716 1.0000000 [110,] 0.043976036 0.043976036 1.0000000 [111,] 0.045529351 0.045529351 1.0000000 [112,] 0.047137531 0.047137531 1.0000000 [113,] 0.048802516 0.048802516 1.0000000 [114,] 0.050526311 0.050526311 1.0000000 [115,] 0.052310993 0.052310993 1.0000000 [116,] 0.054158714 0.054158714 1.0000000 [117,] 0.056071699 0.056071699 1.0000000 [118,] 0.058052255 0.058052255 1.0000000 [119,] 0.060102768 0.060102768 1.0000000 [120,] 0.062225708 0.062225708 1.0000000 [121,] 0.064423635 0.064423635 1.0000000 [122,] 0.066699197 0.066699197 1.0000000 [123,] 0.069055135 0.069055135 1.0000000 [124,] 0.071494290 0.071494290 1.0000000 [125,] 0.074019600 0.074019600 1.0000000 [126,] 0.076634109 0.076634109 1.0000000 [127,] 0.079340967 0.079340967 1.0000000 [128,] 0.082143436 0.082143436 1.0000000 [129,] 0.085044893 0.085044893 1.0000000 [130,] 0.088048836 0.088048836 1.0000000 [131,] 0.091158883 0.091158883 1.0000000 [132,] 0.094378783 0.094378783 1.0000000 [133,] 0.097712415 0.097712415 1.0000000 [134,] 0.101163798 0.101163798 1.0000000 [135,] 0.104737090 0.104737090 1.0000000 [136,] 0.108436597 0.108436597 1.0000000 [137,] 0.112266777 0.112266777 1.0000000 [138,] 0.116232247 0.116232247 1.0000000 [139,] 0.120337784 0.120337784 1.0000000 [140,] 0.124588336 0.124588336 1.0000000 [141,] 0.128989026 0.128989026 1.0000000 [142,] 0.133545156 0.133545156 1.0000000 [143,] 0.138262217 0.138262217 1.0000000 [144,] 0.143145894 0.143145894 1.0000000 [145,] 0.148202071 0.148202071 1.0000000 [146,] 0.153436841 0.153436841 1.0000000 [147,] 0.158856513 0.158856513 1.0000000 [148,] 0.164467618 0.164467618 1.0000000 [149,] 0.170276917 0.170276917 1.0000000 [150,] 0.176291412 0.176291412 1.0000000 [151,] 0.182518349 0.182518349 1.0000000 [152,] 0.188965234 0.188965234 1.0000000 [153,] 0.195639834 0.195639834 1.0000000 [154,] 0.202550194 0.202550194 1.0000000 [155,] 0.209704640 0.209704640 1.0000000 [156,] 0.217111795 0.217111795 1.0000000 [157,] 0.224780583 0.224780583 1.0000000 [158,] 0.232720248 0.232720248 1.0000000 [159,] 0.240940356 0.240940356 1.0000000 [160,] 0.249450814 0.249450814 1.0000000 [161,] 0.258261876 0.258261876 1.0000000 [162,] 0.267384162 0.267384162 1.0000000 [163,] 0.276828663 0.276828663 1.0000000 [164,] 0.286606762 0.286606762 1.0000000 [165,] 0.296730241 0.296730241 1.0000000 [166,] 0.307211300 0.307211300 1.0000000 [167,] 0.318062569 0.318062569 1.0000000 [168,] 0.329297126 0.329297126 1.0000000 [169,] 0.340928507 0.340928507 1.0000000 [170,] 0.352970730 0.352970730 1.0000000 [171,] 0.365438307 0.365438307 1.0000000 [172,] 0.378346262 0.378346262 1.0000000 [173,] 0.391710149 0.391710149 1.0000000 [174,] 0.405546074 0.405546074 1.0000000 [175,] 0.419870708 0.419870708 1.0000000 [176,] 0.434701316 0.434701316 1.0000000 [177,] 0.450055768 0.450055768 1.0000000 [178,] 0.465952567 0.465952568 1.0000000 [179,] 0.482410870 0.482410893 1.0000000 [180,] 0.499450512 0.499450548 1.0000000 [181,] 0.517092024 0.517110629 1.0000000 [182,] 0.535356668 0.534935744 1.0000000 [183,] 0.554266452 1.000000000 1.0000000 [184,] 0.573844165 1.000000000 1.0000000 [185,] 0.594113398 1.000000000 1.0000000 [186,] 0.615098579 1.000000000 1.0000000 [187,] 0.636824994 1.000000000 1.0000000 [188,] 0.659318827 1.000000000 1.0000000 [189,] 0.682607183 1.000000000 1.0000000 [190,] 0.706718127 1.000000000 1.0000000 [191,] 0.731680714 1.000000000 1.0000000 [192,] 0.757525026 1.000000000 1.0000000 [193,] 0.784282206 1.000000000 1.0000000 [194,] 0.811984499 1.000000000 1.0000000 [195,] 0.840665289 1.000000000 1.0000000 [196,] 0.870359136 1.000000000 1.0000000 [197,] 0.901101825 1.000000000 1.0000000 [198,] 0.932930403 1.000000000 1.0000000 [199,] 0.965883224 1.000000000 1.0000000 [200,] 1.000000000 1.000000000 1.0000000 > showProc.time() Time (user system elapsed): 0.17 0.02 0.19 > > > ## Excursion: In order to check if the extreme-left tail approximation is ok, > ## we need to compute log(p*beta(p,q)) accurately for small p > ## This suffers from "cancellation" and is for *EXPERIMENTATION* > ## --> c12pBeta() further below is for "production" > c12.pBeta <- function(q, c2.only=TRUE) { + ## Purpose: Compute 1st and 2nd coefficient of p*Beta(p,q)= 1 + c_1*p + c_2*p^2 + .. + ## ---------------------------------------------------------------------- + ## Arguments: q -- can be a vector + ## ---------------------------------------------------------------------- + ## Author: Martin Maechler, Date: 30 Oct 2009 + stopifnot(q > 0) + psi1 <- digamma(1)# == - (Euler's)gamma = -0.5772157 + psiD1 <- trigamma(1) + + c1 <- psi1 - digamma(q) + ## Note: psi(1)^2 - psi(q)^2 = (psi(1) - psi(q))*(psi(1)+psi(q)) = + ## and (trigamma(1) - trigamma(q) + c1*(psi1+psi(q)))/2 - psi(q) * c1 + ##.... + c2 <- (psiD1 - trigamma(q) + c1^2)/2 + ## Note: cancellation for q -> 0 !! ---> c12pBeta() below + ## digamma(q) = psi (q) ~= -1/q + psi(1) + q*psi'(1) + q^2/2*psi''(1) +... + ## trigamma(q) = psi'(q) ~= 1/q^2 + psi'(1) + q*psi''(1) ... + ## ==> 2*c_2 = psi'(1) - psi'(q) + (psi(q) - psi(1))^2 = + ## = - 1/q^2 -q*psi''(1) -.. + (-1/q + q*psi'(1) + q^2/2*psi''(1)+..)^2 = + ## = -q*psi''(1) -2*psi'(1) - q*psi''(1) = 2*(-psi'(1)-q*psi''(1))+ O(q^2) + ## ==> c_2(q) = -(psi'(1) + q*psi''(1)) + O(q^2) + ## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + if(c2.only) c2 else cbind(c1 = c1, c2 = c2) + } > c12.pBeta(1e-5,FALSE) c1 c2 [1,] 1e+05 -1.644912 > summary(warnings()) 9 identical warnings: In qbetaAppr(p, a, b) : a[.] > 1/2 (not thought for!) > > showProc.time() Time (user system elapsed): 0 0 0 > > ## c1 c2 > ## [1,] 1e+05 -1.644912 > trigamma(1) [1] 1.644934 > ## [1] 1.644934 > ## which confirms the expansion of c_2(q) above > > curve(c12.pBeta, .001, 100) ##-> > curve(c12.pBeta, 1, 1000, log="xy") ## dominated by -digamma(q)^2 / 2 Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 1 y value <= 0 omitted from logarithmic plot > ## Cancellation problem when q -> 0 as well ... probably not really relevant > ## *but* fixed in c12pBeta() below: > curve(c12.pBeta, 5e-9, 100, log="x") > > tit.pBeta <- function() { + title("c12.pBeta(x) and its (log-)linear approximation") + mtext(substitute(c[1] + c[2]*x == C1 + C2*x, + setNames(as.list(cq12), c("C1", "C2"))), col=2) + } > cq12 <- -psigamma(1, 1:2) ## == c( -digamma(1), -trigamma(1) ) > curve(cq12[1] + cq12[2]*x, add=TRUE, col=2, n=1000) > tit.pBeta() > > ## a bit zoomed: > curve(c12.pBeta, 1e-8, .2, log="x") > curve(cq12[1] + cq12[2]*x, add=TRUE, col="red", n=1000); tit.pBeta() > ## difference only: > curve(c12.pBeta(x) - (cq12[1] + cq12[2]*x), + 1e-8, .2, log="x", ylim = .01*c(-1,1)); abline(h=0, lty=3, col="gray60") > ## | . - . | log-scaled: > curve(abs(c12.pBeta(x) - (cq12[1] + cq12[2]*x)), + 1e-8, .2, log="xy") > abline(v= 4e-4, col="lightblue") > ##==> use q = 4e-4 as cutoff > > qq <- lseq(1e-6, 1e-2, length= 64) > 1 - (cq12[1] + cq12[2]*qq) / c12.pBeta(qq) [1] 1.631248e-05 -5.767047e-05 5.391123e-05 -1.435149e-06 -1.963130e-05 [6] -9.941867e-06 9.088572e-06 5.003631e-06 1.005729e-06 6.384822e-06 [11] 2.852923e-07 -1.041129e-06 -1.051746e-06 -8.814963e-07 7.786502e-08 [16] 1.279522e-06 4.425331e-07 7.997872e-07 8.328916e-08 9.150328e-08 [21] 2.321904e-08 -4.233614e-08 -8.593431e-08 1.435097e-07 -6.460389e-08 [26] -2.453752e-08 -6.565086e-09 -5.575289e-08 -2.411003e-08 8.336898e-09 [31] 1.729125e-09 1.253202e-08 1.512903e-08 1.901125e-08 1.969310e-08 [36] 2.591769e-08 2.749001e-08 4.446160e-08 5.411663e-08 7.364086e-08 [41] 9.818004e-08 1.323786e-07 1.769391e-07 2.372900e-07 3.180417e-07 [46] 4.260885e-07 5.710152e-07 7.648694e-07 1.024871e-06 1.373139e-06 [51] 1.839662e-06 2.464802e-06 3.302437e-06 4.424944e-06 5.929124e-06 [56] 7.944924e-06 1.064650e-05 1.426743e-05 1.912092e-05 2.562716e-05 [61] 3.434989e-05 4.604566e-05 6.173006e-05 8.276691e-05 > plot(qq, 1 - c12.pBeta(qq)/(cq12[1] + cq12[2]*qq), type="o", log="x") > plot(qq, abs(1 - c12.pBeta(qq)/(cq12[1] + cq12[2]*qq)), type="o", log="xy") > > showProc.time() Time (user system elapsed): 0.02 0.02 0.03 > > ## Compute log(p*beta(p,q)) accurately for small p > c12pBeta <- function(q, c2.only=FALSE, cutOff.q = 4e-4) { + ## Purpose: Compute 1st and 2nd coefficient of p*Beta(p,q) = 1 + c_1*p + c_2*p^2 + .. + ## NB: The Taylor expansion of log(p*Beta(p,q)) = c_1*p + d_2*p^2 + .. + ## is simpler: c_1=d_1, and d_j = \psi^(j)(1) - \psi^(j)(q) = psigamma(1,j) - psigamma(q,j) + ## ---------------------------------------------------------------------- + ## Arguments: q -- can be a vector + ## ---------------------------------------------------------------------- + ## Author: Martin Maechler, Date: 30 Oct 2009 + stopifnot(q > 0, cutOff.q >= 0) + psi1 <- digamma(1) # == - (Euler's)gamma = -0.5772157 + psiD1 <- trigamma(1)# == pi^2 / 6 + + ## vectorize in q: + r <- q + if(!c2.only) + c1 <- psi1 - digamma(q) + if(any(sml <- q < cutOff.q)) { + ## To avoid cancellation, use Taylor expansion for small q: + ## digamma(q)= psi (q) ~= -1/q + psi(1) + q*psi'(1)+q^2/2*psi''(1)+.. + ## trigamma(q)= psi'(q) ~= 1/q^2+ psi'(1)+ q*psi''(1) ... + ## ==> 2*c_2 = psi'(1) - psi'(q) + (psi(q) - psi(1))^2 = + ## = - 1/q^2 -q*psi''(1) -.. + (-1/q + q*psi'(1) + q^2/2*psi''(1)+..)^2= + ## = -q*psi''(1) -2*psi'(1) - q*psi''(1) = 2*(-psi'(1)-q*psi''(1))+ .. + ## ==> c_2(q) = -(psi'(1) + q*psi''(1)) + O(q^2) + ## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + r[sml] <- -(psiD1 + psigamma(1,2) * q[sml]) + } + if(!all(sml)) { ## regular case: + i.ok <- !sml + q <- q[i.ok] + c1. <- if(c2.only) psi1 - digamma(q) else c1[i.ok] + + ## Note: psi(1)^2 - psi(q)^2 = (psi(1) - psi(q))*(psi(1)+psi(q)) = + ## and (trigamma(1) - trigamma(q) + c1*(psi1+psi(q)))/2 - psi(q) * c1 + ##.... + r[i.ok] <- (psiD1 - trigamma(q) + c1.^2)/2 + } + if(c2.only) r else cbind(c1 = c1, c2 = r) + } > > tit <- quote(c[2](q) ~~~ "from Beta expansion" ~~ p*Beta(p,q) %~~% 1 + c[1]*p + c[2]*p^2) > doX <- function(cutOff.q = eval(formals(c12pBeta)$cutOff.q)) { + eaxis(1); py <- par("usr")[3:4] + abline(h=0, v=1, lty=3) + text(1, py[1]+0.95*diff(py), quote(q==1), adj=-.005) + abline(v = cutOff.q, col=2, lty=2, lwd=2) + text(cutOff.q, py[1]+0.9*diff(py), + paste0("cutOff.q=",formatC(cutOff.q)), col=2, adj=-0.05) + } > curve(c12pBeta(x,TRUE), 5e-9, 100, log="x", n=1025, xaxt='n', + xlab=quote(q), main = tit); doX()#--> no cancellation problems > showProc.time() Time (user system elapsed): 0.03 0 0.03 > > ## zoom in > curve(c12pBeta(x,TRUE), 3e-4, 5e-4, log="x", n=1025, xaxt='n', + xlab=quote(q), main = tit); doX()#--> no break.. > showProc.time() Time (user system elapsed): 0 0 0 > > ## zoom out > curve(c12pBeta(x,TRUE), 1e-100, 1e110, log="x", n=1025, xaxt='n', + xlab=quote(q), main = tit); doX()# .. neither problems .. > > curve(c12pBeta(x,TRUE), 1e-100, 1e110, log="x", n=1025, xaxt='n', ylim = c(-2,1), + xlab=quote(q), main = tit); doX()# .. neither problems .. > > summary(warnings()) 1 identical warnings: In xy.coords(x, y, xlabel, ylabel, log) : 1 y value <= 0 omitted from logarithmic plot > showProc.time() Time (user system elapsed): 0.01 0 0.02 > > ## In order for the first oder approximation p*beta(p,q) = 1 + c_1*p to be ok, > ## that c_2*p^2 << 1, i.e. p^2 < eps.c / c_2 <=> p < sqrt(eps.c / c_2) > ## ===================== > ## --> draw the cutoff line below ( --> 'p.cut') > > p.pBeta <- function(q, p.min=1e-20, p.max=1e-6, n=1000, do.ord2=TRUE, log = "xy", ylim=NULL) + { + ## Purpose: Plot log(p * beta(p,q)) for small p -- to visualize cancellation + ## ---------------------------------------------------------------------- + ## Arguments: + ## ---------------------------------------------------------------------- + ## Author: Martin Maechler, Date: 30 Oct 2009, 09:16 + stopifnot(0 < p.min, p.min < p.max, n > 100, n == round(n), + is.character(log), length(log) == 1, + length(q) == 1, 0 < q) + sgn <- if(q > 1) -1 else 1 + curve(sgn*log(x*beta(x,q)), p.min, p.max, n=n, log=log, ylim=ylim, + col="gray", lwd=3, ylab="", xlab = expression(p), + axes = FALSE, frame.plot = TRUE, # rather use eaxis() below + main = substitute("Accurately computing "~ log(p %.% B(p,q)) ~ + " for "~ {p %->% 0} *", "~ (q == Q), list(Q=format(q,digits=15)))) + eaxis(1); eaxis(2) + curve(sgn*(log(x)+lbeta(x,q)), add = TRUE, col="red3", n=n) + ## Here comes the 1st order cancellation solution + c1 <- sgn*( digamma(1) - digamma(q)) + c2 <- sgn*(trigamma(1) - trigamma(q))/2 + colc1 <- adjustcolor("skyblue", 1/2) + c1y <- curve(c1 * x, add = TRUE, col=colc1, lwd=2, n=n)$y + ## the 2nd order approx. + if(do.ord2) { + colc2 <- "#60C07064" # opaque (!) + c12y <- curve(x*(c1 + c2* x), add = TRUE, col=colc2, lwd=3, n=n)$y + } + ## rather: + collog1 <- adjustcolor("forestgreen", 1/2) + lc1y <- curve(sgn*log1p(sgn*c1*x), add = TRUE, col=collog1, lty=5, lwd=3, n=n)$y + ## the cutoff, from where 1st order should be perfect: + p.cut <- sqrt(.Machine$double.eps / abs(c2)) + txt.cut <- "cutoff(1st order approx.)" + cat(txt.cut,": ", formatC(p.cut),"\n", sep="") + ## FIXME: These *badly* fail, if we do NOT have log = "xy" : + ux <- par("xaxp")[1:2] # e.g. [1e-20, 1e-6] + uy <- par("yaxp")[1:2] + in.x <- function(f) { stopifnot(0 <= f, f <= 1); 10^c(c(1-f, f) %*% log10(ux)) } + in.y <- function(f) { stopifnot(0 <= f, f <= 1); 10^c(c(1-f, f) %*% log10(uy)) } + if(ux[1] <= p.cut & p.cut <= ux[2]) { + abline(v = p.cut, col = "blue", lwd=2, lty=2) + text(p.cut, in.y(.96), txt.cut, col="blue") + } + else if(p.cut > ux[2]) text(in.x(.98), in.y(0.5), paste("-->", txt.cut), col = "blue", lwd=2, lty=2) + else if(p.cut < ux[1]) text(in.x(.02), in.y(0.5), paste(txt.cut, "<--"), col = "blue", lwd=2, lty=2) + leg <- + if(sgn == -1) expression(-log(p*B(p,q)), -(log(p) + log(B(p,q))), -(psi(1)-psi(q))*p, -log1p(c[1]*p)) + else expression(log(p*B(p,q)), log(p) + log(B(p,q)), (psi(1)-psi(q))*p, log1p(c[1](q)*p)) + if(do.ord2) { leg <- leg[c(1:4,4)]; leg[[4]] <- quote(p*(c[1](q) + c[2](q)* p)) } + legend("topleft", legend= leg, inset=.01, + lwd= c(3,1,2,if(do.ord2)3, 3), + col= c("gray","red3", colc1, if(do.ord2)colc2, collog1)) + mtext(R.version.string, side=4, cex = 3/5, adj=0) + } > > showProc.time() Time (user system elapsed): 0.02 0 0.01 > > p.pBeta(q = 1.1) ##q should not matter much; q very close to 1 will be another challenge cutoff(1st order approx.): 4.581e-08 Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 132 y values <= 0 omitted from logarithmic plot > ## zoom in: > p.pBeta(1.1, p.max = 1e-12) cutoff(1st order approx.): 4.581e-08 Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 263 y values <= 0 omitted from logarithmic plot > p.pBeta( 11, p.max = 1e-11)## log(..) is not much worse here cutoff(1st order approx.): 1.693e-08 Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 157 y values <= 0 omitted from logarithmic plot > p.pBeta(101, p.max = 1e-11)## !! interestingly, here log(..) is not worse cutoff(1st order approx.): 1.648e-08 Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 308 y values <= 0 omitted from logarithmic plot > p.pBeta(1001, p.max = 1e-11)## both versions are identical cutoff(1st order approx.): 1.644e-08 Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 393 y values <= 0 omitted from logarithmic plot > > showProc.time() Time (user system elapsed): 0.11 0 0.11 > ## q -> 1 > p.pBeta(1.001, p.max = 1e-10) cutoff(1st order approx.): 4.301e-07 Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 540 y values <= 0 omitted from logarithmic plot > p.pBeta(1+1e-5, 1e-15, 1e-8)# problems start earlier! cutoff(1st order approx.): 4.298e-06 Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 225 y values <= 0 omitted from logarithmic plot > p.pBeta(1+1e-7, 1e-11, 1e-6)# problems start even earlier! cutoff(1st order approx.): 4.298e-05 Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 252 y values <= 0 omitted from logarithmic plot > p.pBeta(1+1e-9, 1e-9, 1e-4)# problems start even earlier! cutoff(1st order approx.): 0.0004298 Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 257 y values <= 0 omitted from logarithmic plot > ## > showProc.time() Time (user system elapsed): 0.05 0 0.05 > ## q < 1 > p.pBeta(0.9) cutoff(1st order approx.): 4e-08 Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 357 y values <= 0 omitted from logarithmic plot > p.pBeta(0.1, p.max = 1e-8) cutoff(1st order approx.): 2.11e-09 Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 277 y values <= 0 omitted from logarithmic plot > p.pBeta(0.01) cutoff(1st order approx.): 2.107e-10 Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 140 y values <= 0 omitted from logarithmic plot > p.pBeta(0.0001) cutoff(1st order approx.): 2.107e-12 Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 4 y values <= 0 omitted from logarithmic plot > p.pBeta(0.0001,1e-23, 1e-4) cutoff(1st order approx.): 2.107e-12 Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 150 y values <= 0 omitted from logarithmic plot > p.pBeta(0.0001,1e-12, 1e-2) cutoff(1st order approx.): 2.107e-12 > p.pBeta(0.0001,1e-6, 1e-3) cutoff(1st order approx.): 2.107e-12 > p.pBeta(0.0001,1e-6, 1e-1) cutoff(1st order approx.): 2.107e-12 > > summary(warnings()) 1 identical warnings: In xy.coords(x, y, xlabel, ylabel, log) : 150 y values <= 0 omitted from logarithmic plot > showProc.time() Time (user system elapsed): 0.09 0.01 0.11 > > if(!require("Rmpfr")) quit("no") Loading required package: Rmpfr Loading required package: gmp Attaching package: 'gmp' The following objects are masked from 'package:sfsmisc': factorize, is.whole The following objects are masked from 'package:base': %*%, apply, crossprod, matrix, tcrossprod C code of R package 'Rmpfr': GMP using 64 bits per limb Attaching package: 'Rmpfr' The following object is masked from 'package:gmp': outer The following object is masked from 'package:DPQ': log1mexp The following objects are masked from 'package:stats': dbinom, dgamma, dnbinom, dnorm, dpois, dt, pnorm The following objects are masked from 'package:base': cbind, pmax, pmin, rbind > ##--===============---====------------ Rmpfr needed from here ----------------------------- > > > p.err.pBeta <- function(q, p.min=1e-12, p.max=1e-6, n=1000, + kind = c("relErr", "absErr", "error")) + { + ## Purpose: Plot log(p * beta(p,q)) for small p -- to visualize cancellation + ## ---------------------------------------------------------------------- + ## Arguments: + ## ---------------------------------------------------------------------- + ## Author: Martin Maechler, Date: 30 Oct 2009, 09:16 + stopifnot(0 < p.min, p.min < p.max, n > 100, n == round(n), + length(q) == 1, 0 < q) + + c12 <- c12pBeta(q) + d2 <- (trigamma(1) - trigamma(q))/2 + c1 <- c12[1] ## == digamma(1) - digamma(q) + c2 <- c12[2] + + p <- lseq(p.min, p.max, length=n) + + T0 <- log(p*beta(p,q))# which we know has its problems, too + T1 <- c1*p + T2 <- p*(c1 + d2*p) + T3 <- log1p(T1) + T4 <- log1p(p*(c1 + c2*p)) + + ## Better: compute "true value" very accurately: + stopifnot(require("Rmpfr")) + p. <- mpfr(p, prec=200) + tt <- log(p.*beta(p.,q)) + + p.cut <- sqrt(.Machine$double.eps / abs(c2)) + cat("cutoff :", formatC(p.cut),"\n") + + err.mat <- as(tt - cbind(T0, T1, T3, T2, T4), + "matrix")## classic numeric matrix + leg <- expression(T - log(p*beta(p,q)) ~ " (double)", + T - c[1]*p, # T1 + T - log1p(c[1]*p), # T3 + T - (c[1]*p+d[2]*p^2), # T2 + T - log1p(c[1]*p+c[2]*p^2))# T4 + stopifnot(length(leg) == ncol(err.mat)) + main <- paste("log(p*beta(p,q)) approx. errors, q = ", formatC(q)) + ## `col` below is from + ## dput(RColorBrewer::brewer.pal(length(leg), "Dark2")) + col <- c("#1B9E77", "#D95F02", "#E7298A", "#7570B3", "#66A61E") + ## put the color which is "most gray" (i.e equal R,G,B values) first: + nc <- length(col) + cm <- col2rgb(col) + i1 <- which.min(colSums(sweep(cm, 2, colMeans(cm))^2)) + if(i1 > 1) col <- col[c(i1:nc, 1:(i1-1))] + m.matplot <- function(x,y, log, main, ..., LEG, xy.leg) { + log.. <- strsplit(log, "")[[1]] + log.x <- any("x" == log..) + log.y <- any("y" == log..) + + matplot(x,y, log=log, ..., xlab = quote(p), type = "l", col=col, main=main, + xaxt = if(log.x) "n" else "s", + yaxt = if(log.y) "n" else "s") + if(log.x) eaxis(1, nintLog = 15) + if(log.y) eaxis(2, nintLog = 20) # "FIXME" eaxis(): default nintLog=10 does not cover well + legend(xy.leg, LEG, lty=1:5, col=col, inset=.01) + } + kind <- match.arg(kind) + switch(kind, + "error" = + m.matplot(p, err.mat, log="x", main=main, + ylab = "log(p*beta(p,q)) - 'approx.'", + LEG = leg, xy.leg = "left"), + "absErr" = { + + m.matplot(p, abs(err.mat), ## <- absolute -- log scale + log="xy", main=main, + ylab = "|log(p*beta(p,q)) - 'approx.'|", + LEG = as.expression(lapply(leg, + function(.) substitute(abs(EXPR), list(EXPR= .)))), + xy.leg = "topleft") + }, + "relErr" = { + m.matplot(p, abs(err.mat/as(tt,"numeric")), ## <- RELATIVE errors + log="xy", main = paste( + "log(p*beta(p,q)) RELATIVE approx. errors, q = ", + formatC(q)), + ylab = "|1 - 'approx'/log(p*beta(p,q))|", + LEG = as.expression(lapply(leg, + function(.) substitute(abs(EXPR), list(EXPR= .)))), + xy.leg = "topleft") + }, + stop("invalid 'kind' : ", kind) + )# end{switch} + + abline(v = p.cut, col = "blue", lwd=2, lty=2) + mtext(R.version.string, side=4, cex = 3/5, adj=0) + } > > ## with q >> 1, the log1p() approx. are *worse* (see more below) ! > p.err.pBeta(21, 1e-18) # cutoff : 5.527e-09 cutoff : 5.527e-09 > p.err.pBeta(21, , 1e-3) cutoff : 5.527e-09 > summary(warnings()) 1 identical warnings: In xy.coords(x, y, xlabel, ylabel, log) : 150 y values <= 0 omitted from logarithmic plot > showProc.time() Time (user system elapsed): 0.64 0.03 0.67 > > if(doExtras) { # each plot is somewhat large & expensive ..------------------- + ## hmm, with q >> 1, the log1p() approx. are *worse* .. + p.err.pBeta(q = 2.1, kind="absErr") + p.err.pBeta(q = 2.1, 1e-13, kind="absErr") + p.err.pBeta(1.1, 1e-13, 0.1)## log1p() still very slightly *worse* + ## here, q < 1 and log1p() are already (slightly) better + p.err.pBeta(0.8) + p.err.pBeta(0.8, 1e-10, 1e-2) + p.err.pBeta(0.1) # log1p(<2 trms>) is best; zoom out (to guess "cutoff_2") + p.err.pBeta(0.1, 1e-15, 1e-3) + p.err.pBeta(0.001,, 1e-2)## Here, log1p(<2 terms>) is clearly best + ## reminding of the test qbeta(1e-300, 2^-12, 2^-10): + p.err.pBeta(2^-10, p.max = 2^-10); abline(v=2^-12, lty=3, col="gray20") + ## --> ok, the (*, 2^-12, 2^-10) is *not* yet critical + + ##---> q < 1 ==> the log1p() versions are superior + ##---> q > 1 ==> the direct versions are superior + ## in both cases: The quadratic term is an order of magnitude better + ## the above p.cut is "order-of-magnitude" correct, but not really numerically ok! + + print(summary(warnings())) + showProc.time() + }# only if(doExtras) ---------------------------------------------------------------------- > > pBeta <- function(p, q) + { + ## Purpose: Compute log(p * beta(p,q)) accurately also for small p + ## ---------------------------------------------------------------------- + ## Arguments: p: (vector) q: scalar + ## ---------------------------------------------------------------------- + ## Author: Martin Maechler, Date: 31 Oct 2009, 18:38 + + stopifnot(0 < p, 0 < q, length(q) == 1, (n <- length(p)) > 0) + + if(q == 1) return(0*p) + + c12 <- c12pBeta(q) + d2 <- (trigamma(1) - trigamma(q))/2 + c1 <- c12[1] ## == digamma(1) - digamma(q) + c2 <- c12[2] + + ## the cutoff from where the approximation is "perfect" + p.cut <- sqrt(.Machine$double.eps / abs(d2)) + r <- p + if(any(sml <- p <= p.cut)) { + p. <- p[sml] + r[sml] <- if(q < 1) log1p(p.*(c1 + c2*p.)) else p.*(c1 + d2*p.) + } + if(any(ok <- !sml)) { + p <- p[ok] + r[ok] <- log(p * beta(p, q)) + } + r + } > > ### More on log(p * Beta(p,q)) for p << 1 > ### ================== > > > ## 1. Slightly nicer Gamma(.) plot than in example(Gamma) : ---------- > op <- options(warn = -1) > x <- sort(c(seq(-3, 4, length.out = 201), outer(0:-3, (-1:1)*1e-6, "+"))) > plot(x, gamma(x), ylim = c(-20,20), col = "red", type = "l", lwd = 2, + main = expression(Gamma(x))) > abline(h = 0, v = -3:0, lty = 3, col = "midnightblue") > ## Warning .. in gamma(x) : NaNs produced > axis(4, at=factorial(1:3), las=2, cex.axis=3/4) > abline(h=1, lty=4, col="thistle") > for(t in c("p","h")) points(1:4, gamma(1:4), pch=4, type=t, lty=2) > > ## Now try, comparing with Rmpfr exact values: > if(!require("Rmpfr")) { message("Cannot attach 'Rmpfr' package"); q("no") } > > a. <- 1/256 # is exact! > b. <- 1/2 # " " > a <- mpfr(a., 512) > b <- mpfr(b., 512) > (aBa1 <- log(a*beta(a,b))) 1 'mpfr' number of precision 512 bits [1] 0.00539025506615457152696124105278305377723163029605479125314934854012553628647620193253637462441295398525816601026739203604439887171688979382925584525277609414 > ## 0.0053902550661545715269612410527830537772316302960547912531493485401255362864762019.... > (aBa2 <- log(a)+lbeta(a,b))# 1 'mpfr' number of precision 512 bits [1] 0.00539025506615457152696124105278305377723163029605479125314934854012553628647620193253637462441295398525816601026739203604439887171688979382925584525277627886 > ## 0.0053902550661545715269612410527830537772316302960547912531493485401255362864762019.... > ## how close? > asNumeric((aBa1-aBa2)*2/(aBa1+aBa2)) # -3.426748e-152; (NB 512 bits ~= 154.1 digits) [1] -3.426748e-152 > > ## Double precison accuracy > dput(log(a. * beta(a.,b.))) 0.00539025506615481 > ## 0.00539025506615481 > dput(log(a.)+lbeta(a.,b.))# (to my surprise!) slightly worse 0.00539025506615509 > ## 0.00539025506615509 > asNumeric( log(a. * beta(a.,b.))/(log(a)+lbeta(a,b)) - 1) # 4.504144e-14 [1] 4.504144e-14 > asNumeric((log(a.)+lbeta(a.,b.))/(log(a)+lbeta(a,b)) - 1) # 9.653358e-14 [1] 9.653358e-14 > > psigamma(1) # -0.5772157 = Euler's gamma [1] -0.5772157 > > ## Use my approx: > a.* (psigamma(1) - psigamma(b.)) ## 0.005415212 [1] 0.005415212 > a.* (psigamma(1) - psigamma(b.) + a./2*(psigamma(1, d=1) - psigamma(b., d=1))) [1] 0.005390113 > ## 0.005390113 {clearly better, still only ~ 4 correct digits > > showProc.time() Time (user system elapsed): 0.09 0 0.1 > > proc.time() user system elapsed 2.71 0.29 3.00