R Under development (unstable) (2024-08-17 r87027 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. > require(DPQmpfr) Loading required package: DPQmpfr > require("Rmpfr") # (is in DPQmpfr's strict dependencies) Loading required package: Rmpfr Loading required package: gmp Attaching package: 'gmp' The following objects are masked from 'package:base': %*%, apply, crossprod, matrix, tcrossprod C code of R package 'Rmpfr': GMP using 64 bits per limb Attaching package: 'Rmpfr' The following object is masked from 'package:gmp': outer The following objects are masked from 'package:stats': dbinom, dgamma, dnbinom, dnorm, dpois, dt, pnorm The following objects are masked from 'package:base': cbind, pmax, pmin, rbind > > options(warn = 1)# warnings *immediately* > > ## *.time() utilities from > source(system.file(package="Matrix", "test-tools-1.R", mustWork=TRUE)) Loading required package: tools > source(system.file(package="DPQ", "test-tools.R", mustWork=TRUE)) > ## => list_() , loadList() , readRDS_() , save2RDS() > > (doExtras <- DPQmpfr:::doExtras()) [1] FALSE > ## save directory (to read from): > (sdir <- system.file("safe", package="DPQmpfr")) [1] "D:/RCompile/CRANincoming/R-devel/lib/DPQmpfr/safe" > ## initially, when we have old version of pkg (only for pkg maintainer!!) > if(!nzchar(sdir) && dir.exists(pDir <- "~/R/Pkgs/DPQmpfr")) { + sdir <- file.path(pDir, "inst/safe") + if(!dir.exists(sdir)) dir.create(sdir) + } > ## IGNORE_RDIFF_BEGIN > (has.sdir <- dir.exists(sdir)) [1] TRUE > if(!has.sdir) + cat("safe directory ", sQuote(sdir), " does not exist!\n") > ## IGNORE_RDIFF_END > > > showProc.time(,1) # 1: only 'user' Time (user): 0.02 > > ### pbeta*() : ------------------------------------------------------------- > show_pbD94 <- function(...) showSys.time(show(pbetaD94(..., verbose = TRUE))) > ## ~~~~~~~~ > > show_pbD94(mpfr(0.864, 140), 5, 5, 54, eps=1e-10)# n=233 at convergence, 0.456302619295.... log_scale = FALSE n=27; F=0.056554259575 after A3-A4 loop n= 233 at convergence 1 'mpfr' number of precision 140 bits [1] 0.45630261929567235741833106018017381811792346 Time user system elapsed Time 1.20 0.00 1.21 > show_pbD94(mpfr(0.864, 140), 5, 5, 54, eps=1e-30)# n=572 " " , 0.45630261933697895485.. log_scale = FALSE n=27; F=0.05655425957495302040101592569982 after A3-A4 loop n= 572 at convergence 1 'mpfr' number of precision 140 bits [1] 0.45630261933697895485422846841559025615080321 Time user system elapsed Time 3.47 0.00 3.47 > > show_pbD94(mpfr(0.956, 140), 5, 5, 170, eps=1e-10)# n= 765, 0.602242164945 log_scale = FALSE n=104; F=0.16132149057 after A3-A4 loop n= 765 at convergence 1 'mpfr' number of precision 140 bits [1] 0.60224216494518105841370442668188481299701479 Time user system elapsed Time 4.41 0.10 4.50 > show_pbD94(mpfr(0.956, 140), 5, 5, 170, eps=1e-20)# n=1323, 0.60224216500116548 log_scale = FALSE n=104; F=0.161321490569710714068 after A3-A4 loop n= 1323 at convergence 1 'mpfr' number of precision 140 bits [1] 0.60224216500116548068143220795401886595112736 Time user system elapsed Time 8.17 0.01 8.19 > > show_pbD94(mpfr(0.8686, 140), 10,10, 54, eps=1e-10)# n=304, 0.918779110843.. log_scale = FALSE n=57; F=0.46793660527 after A3-A4 loop n= 304 at convergence 1 'mpfr' number of precision 140 bits [1] 0.91877911084301486831470044375909570846727322 Time user system elapsed Time 1.25 0.02 1.26 > show_pbD94(mpfr(0.8686, 140), 10,10, 54, eps=1e-20)# n=497, 0.9187791109260769059699 log_scale = FALSE n=57; F=0.467936605273279879678 after A3-A4 loop n= 497 at convergence 1 'mpfr' number of precision 140 bits [1] 0.91877911092607690596990611240335653542050305 Time user system elapsed Time 3.06 0.03 3.09 > > show_pbD94(mpfr(0.8787, 140), 20,20, 54, eps=1e-10)# n=456, 0.99986765729631163 log_scale = FALSE n=125; F=0.52529668439 after A3-A4 loop n= 456 at convergence 1 'mpfr' number of precision 140 bits [1] 0.99986765729631163361241717025738921815379759 Time user system elapsed Time 2.68 0.04 2.71 > show_pbD94(mpfr(0.8787, 180), 20,20, 54, eps=1e-20)# n= log_scale = FALSE n=125; F=0.525296684388920543995 after A3-A4 loop n= 691 at convergence 1 'mpfr' number of precision 180 bits [1] 0.99986765738881454716792652583554690089134775739101043411 Time user system elapsed Time 3.85 0.04 3.88 > > show_pbD94(mpfr(0.922, 140), 20,20, 250, eps=1e-10)# n=743, 0.964119072835964766 log_scale = FALSE n=217; F=0.49453493179 after A3-A4 loop n= 744 at convergence 1 'mpfr' number of precision 140 bits [1] 0.96411907284123334443549412757175046574452552 Time user system elapsed Time 4.49 0.03 4.52 > show_pbD94(mpfr(0.922, 180), 20,20, 250, eps=1e-20)# n= log_scale = FALSE n=217; F=0.494534931791177172872 after A3-A4 loop n= 1118 at convergence 1 'mpfr' number of precision 180 bits [1] 0.96411907293079995228995440485426653017583863136629139638 Time user system elapsed Time 6.61 0.08 6.71 > > if(doExtras) withAutoprint({ + ## MM: Try large lambdas ==> what we need is a *RELATIVE* bound --> now implemented + show_pbD94(mpfr(0.93, 160), 20,25, 10000) # n= 5175 1.008067955986..e-116 + show_pbD94(mpfr(0.93, 160), 20,25, 10000, eps=1e-20) # n= 5513 1.008067956082285281372...e-116 + + show_pbD94(mpfr(0.93, 256), 20,25, 10000, eps=1e-30) # n= 5850 1.008067956082285281381936816362733357442687434229...e-116 + ## double the precBits only + show_pbD94(mpfr(0.93, 512), 20,25, 10000, eps=1e-30) # n= 5850 1.008067956082285281381936816362733357442687434229...e-116 + ## double the precBits again + eps !! + show_pbD94(mpfr(0.93, 1024), 20,25, 10000, eps=1e-50) # n= 6520 1.008067956082285281381936816363684457888595749378086614..e-116 + ## and much more extreme accuracy.. seems to work, too: + show_pbD94(mpfr(0.93, 1024), 20,25, 10000, eps=1e-150)# n= 9828 1.00806795608228528138193681636368445788859574937809624269896.......e-116 + }) > > showProc.time() Time (user system elapsed): 39.73 0.35 40.1 > > ### dbeta*() : > > " TODO " [1] " TODO " > > ### qbeta*() --- careful about speed / time usage > > ## This now works, thanks to log_scale ! > (q.5.250 <- qbetaD94(0.95, .5, 250, ncp = 0)) [1] 0.007661102 attr(,"conv") [1] TRUE attr(,"iterNewton") [1] 9 attr(,"n.inner") [1] 433 171 91 55 36 26 19 19 19 > > ## but using mpfr of course works: > showProc.time() Time (user system elapsed): 0 0.01 0.01 > qbetaD94(1 - 1/mpfr(20 ,64), .5, 250) # 0.00766110246787361539366 1 'mpfr' number of precision 64 bits [1] 0.00766110246787362147663 > qbetaD94(1 - 1/mpfr(20,128), .5, 250, eps = 1e-30, delta=1e-20) # 0.0076611024678614272722717848... 1 'mpfr' number of precision 128 bits [1] 0.007661102467861427272271784849421972046952 > showProc.time() Time (user system elapsed): 17.42 0.14 17.56 > > > qb_sfile <- function(p, precB, mpfr) { + stopifnot(is.numeric(precB), precB == trunc(precB), precB >= 10, + 0 <= p, p < 1, is.logical(mpfr)) + paste0("tests_qbD94_", as.character(precB),"_p", + as.character(round(100*asNumeric(p))), + if(mpfr)"_mpfr", ".rds") + } > > > ## Compare with Table 3 of Baharev_et_al 2017 %% ===> ../man/qbBaha2017.Rd <<<<<<<<<<<< > qbetaD94sim <- function(p = 0.95, # p = 1 - alpha + delta = NULL, + precBits = 128, + saveDir, + aa = c(0.5, 1, 1.5, 2, 2.5, 3, 5, 10, 25), + bb = c(1:15, 10*c(2:5, 10, 25, 50))) + { + stopifnot(dir.exists(saveDir)) + r <- lapply(c(double=FALSE, mpfr=TRUE), function(simMpfr) { ## do *both*: without / with mpfr(.) + cat("simMpfr =", simMpfr, ":\n") + sFile <- file.path(saveDir, qb_sfile(p, precBits, simMpfr)) + if(file.exists(sFile) && (simMpfr || !doExtras)) { + ssR_l <- readRDS_(sFile) + str(ssR_l) + loadList(ssR_l, envir = environment()) + cat("[Ok]\n") + } else { ## do run the simulation -- always if(!simMpfr & doExtras) : + qbet <- matrix(NA_real_, length(aa), length(bb), + dimnames = list(a = formatC(aa), b = formatC(bb))) + if(!simMpfr) { + ## p <- 0.95 # = 1 - alpha <<<--- using double precision everywhere below + ## delta <- 1e-7 # <=> eps = delta^2 = 1e-14 + } else { ## simMpfr: takes *VERY LONG* (~ 60 minutes!) in the loop below ((why ???)) + ## p <- 1 - 1/mpfr(20,128) # = 1 - alpha <<<--- using MPFR everywhere below + ## delta <- 1e-18 + p <- mpfr(p, precBits=precBits) + qbet <- as(qbet, "mpfr") + } + if(is.null(delta)) + delta <- if(simMpfr) 1e-18 else 1e-7 + else stopifnot(is.numeric(delta), delta >= 0) + eps <- delta^2 # default; (delta,eps) are needed for safe-file name + showSys.time( + for(ia in seq_along(aa)) { + a <- aa[ia]; cat("\na=",a,": b=") + for(ib in seq_along(bb)) { + b <- bb[ib]; cat(b," ") + qbet[ia, ib] <- qbetaD94(p, a, b, ncp = 0, delta=delta)# eps := delta^2 + } ##~~~~ ~~~~~~~~ + cat("\n") + } + ) + print(summary(warnings())) # many warnings from qbeta() inaccuracies + ## save2RDS() writes .. + save2RDS(list_(aa, bb, qbet), file = sFile) + cat("[Ok]\n") + } ## --- do run simulation ---------------------------- 0 + qbet + }) # F (simMpfr) + ## return + c(list_(aa, bb), r) + } > > if(interactive()) ## try it out small + rdummy <- qbetaD94sim(a=1, b=1:2, saveDir = tempdir()) > > ## IGNORE_RDIFF_BEGIN > if(has.sdir || doExtras) + r.95 <- qbetaD94sim(p = 0.95, saveDir = if(has.sdir) sdir else tempdir()) simMpfr = FALSE : Reading from D:/RCompile/CRANincoming/R-devel/lib/DPQmpfr/safe/tests_qbD94_128_p95.rds Time (user system elapsed): 0.02 0 0.02 List of 3 $ aa : num [1:9] 0.5 1 1.5 2 2.5 3 5 10 25 $ bb : num [1:22] 1 2 3 4 5 6 7 8 9 10 ... $ qbet: num [1:9, 1:22] 0.903 0.95 0.966 0.975 0.98 ... ..- attr(*, "dimnames")=List of 2 .. ..$ a: chr [1:9] "0.5" "1" "1.5" "2" ... .. ..$ b: chr [1:22] "1" "2" "3" "4" ... [Ok] simMpfr = TRUE : Reading from D:/RCompile/CRANincoming/R-devel/lib/DPQmpfr/safe/tests_qbD94_128_p95_mpfr.rds Time (user system elapsed): 0.01 0 0.02 List of 3 $ aa : num [1:9] 0.5 1 1.5 2 2.5 3 5 10 25 $ bb : num [1:22] 1 2 3 4 5 6 7 8 9 10 ... $ qbet:Class 'mpfrMatrix' [package "Rmpfr"] of dimension c(9, 22) and precision 128 .. 0.902 0.95 0.966 0.975 ... [Ok] > ## IGNORE_RDIFF_END > str(r.95) List of 4 $ aa : num [1:9] 0.5 1 1.5 2 2.5 3 5 10 25 $ bb : num [1:22] 1 2 3 4 5 6 7 8 9 10 ... $ double: num [1:9, 1:22] 0.903 0.95 0.966 0.975 0.98 ... ..- attr(*, "dimnames")=List of 2 .. ..$ a: chr [1:9] "0.5" "1" "1.5" "2" ... .. ..$ b: chr [1:22] "1" "2" "3" "4" ... $ mpfr :Class 'mpfrMatrix' [package "Rmpfr"] of dimension c(9, 22) and precision 128 .. 0.902 0.95 0.966 0.975 ... > ## $ aa > ## $ bb > ## $ double [qbet] > ## $ mpfr [qbet] > loadList(r.95) > rm(double, mpfr) # they shadow R functions > > ## number of correct digits [mpfr is only slightly better -- ??} > data(qbBaha2017, package="DPQmpfr") > dm <- with(r.95, c(length(aa), length(bb))) > stopifnot(exprs = { + is.matrix(qbBaha2017) + identical(dm, dim(qbBaha2017)) + identical(dm, dim(r.95$double)) + identical(dm, dim(r.95$ mpfr )) + }) > > ver <- setNames(, names(r.95)[3:4]) > nD <- lapply(ver, function(nm) asNumeric(-log10(abs(1 - r.95[[nm]] / qbBaha2017)))) > > lapply(nD, round, digits=1) $double b a 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 20 30 40 0.5 13.8 6.6 6.0 6.1 5.9 5.9 6.0 5.8 6.6 5.8 5.7 5.6 5.6 7.0 5.5 6.5 6.1 6.5 1 14.1 6.6 6.6 6.4 6.2 6.2 6.0 7.2 6.2 5.8 6.1 6.1 5.8 6.0 5.8 5.6 6.3 6.5 1.5 6.3 8.2 6.2 6.5 6.1 6.1 6.1 6.5 5.9 6.0 5.8 6.0 5.8 5.9 5.7 5.8 5.8 6.4 2 6.4 6.4 6.3 6.3 6.2 7.3 6.1 6.0 6.1 5.9 6.0 6.1 6.1 6.2 5.8 5.7 5.6 5.9 2.5 6.6 6.6 7.1 6.2 8.5 8.8 6.9 6.8 6.6 6.5 7.2 6.2 6.0 7.1 6.1 5.9 5.6 6.1 3 6.4 6.3 6.3 6.3 6.2 6.2 6.2 6.2 6.4 6.0 6.0 6.1 6.2 6.2 6.7 5.7 5.6 5.5 5 6.7 6.9 6.3 7.2 6.3 6.5 6.1 6.3 6.2 6.1 6.7 6.1 6.1 6.3 6.3 6.0 6.2 5.6 10 6.7 6.6 6.6 6.6 6.4 6.4 6.3 6.2 7.0 6.2 6.2 6.7 6.4 6.5 6.9 6.7 6.6 5.8 25 6.4 7.3 6.4 6.7 7.1 7.2 6.6 7.1 6.5 6.4 6.5 6.7 6.7 6.6 6.3 6.6 6.1 6.9 b a 50 100 250 500 0.5 6.3 6.7 6.5 6.1 1 6.4 5.8 5.8 6.6 1.5 6.3 7.3 6.2 6.7 2 6.4 6.0 6.2 6.4 2.5 5.8 6.2 8.0 5.6 3 5.4 6.2 6.2 5.6 5 5.6 6.4 7.5 5.6 10 6.0 5.6 6.2 6.5 25 6.2 6.2 5.7 6.9 $mpfr b a 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 20 30 40 0.5 16.2 6.6 6.0 6.1 5.9 5.9 6.0 5.8 6.6 5.8 5.7 5.6 5.6 7.0 5.5 6.5 6.1 6.5 1 36.0 6.6 6.6 6.4 6.2 6.2 6.0 7.2 6.2 5.8 6.1 6.1 5.8 6.0 5.8 5.6 6.3 6.5 1.5 6.3 8.2 6.2 6.5 6.1 6.1 6.1 6.5 5.9 6.0 5.8 6.0 5.8 5.9 5.7 5.8 5.8 6.4 2 6.4 6.4 6.3 6.3 6.2 7.3 6.1 6.0 6.1 5.9 6.0 6.1 6.1 6.2 5.8 5.7 5.6 5.9 2.5 6.6 6.6 7.1 6.2 8.5 8.8 6.9 6.8 6.6 6.5 7.2 6.2 6.0 7.1 6.1 5.9 5.6 6.1 3 6.4 6.3 6.3 6.3 6.2 6.2 6.2 6.2 6.4 6.0 6.0 6.1 6.2 6.2 6.7 5.7 5.6 5.5 5 6.7 6.9 6.3 7.2 6.3 6.5 6.1 6.3 6.2 6.1 6.7 6.1 6.1 6.3 6.3 6.0 6.2 5.6 10 6.7 6.6 6.6 6.6 6.4 6.4 6.3 6.2 7.0 6.2 6.2 6.7 6.4 6.5 6.9 6.7 6.6 5.8 25 6.4 7.3 6.4 6.7 7.1 7.2 6.6 7.1 6.5 6.4 6.5 6.7 6.7 6.6 6.3 6.6 6.1 6.9 b a 50 100 250 500 0.5 6.3 6.7 6.5 6.1 1 6.4 5.8 5.8 6.6 1.5 6.3 7.3 6.2 6.7 2 6.4 6.0 6.2 6.4 2.5 5.8 6.2 8.0 5.6 3 5.4 6.2 6.2 5.6 5 5.6 6.4 7.5 5.6 10 6.0 5.6 6.2 6.5 25 6.2 6.2 5.7 6.9 > ## well, Ding(1994) may not be so good? > > matplot (bb, t(nD$ double), type = "o", log="xy", xaxt="n") > matlines(bb, t(nD$ mpfr )) # again: mpfr is almost the same > axis(1, at=bb, cex.axis = 0.8, mgp = (2:0)/2) > > showProc.time() Time (user system elapsed): 0.11 0 0.1 > > ## Looks not good: > > ## NB: *relative* convergence is not good here for f() ~= 0 !!! <<< Ding *did* have absolute > ## ==> use same idea as 'nls(... control = list(scaleOffset=1))' !! > > a <- 1.5 > b <- 7 > > if(interactive())## verbose=3 now suggest it's hopeless (?!!) + try({ ## the 'mpfr' version *does* converge "n = 372" + qbb <- list(mpfr = qbetaD94(mpfr(0.95, 128), shape1= a, shape2 = b, ncp = 100, delta=1e-18, itrmax = 2000, verbose=3), + dble = qbetaD94(0.95, shape1= a, shape2 = b, ncp = 100, delta=1e-7, itrmax = 600, verbose=3)) + ## also itrmax = 1e6 fails the same + ## even delta = 1e-2 fails: reason: f() = F'() = 1.25e-59 <<< bound2 ~= 3e-16 + }) > > ## Note how silly it is to have a very small 'eps' in a situation were 'x' is still far from the truth > ## ===> Idea: Much faster if 'eps' is "large" at the beginning, when the Newton 'd.x' will be inaccurate anyway !! > > ## Before we had 'log_scale', ncp=100 could not work > pbetaD94(.99, a, b, ncp = 100, log_scale=FALSE, verbose=TRUE) -> pb. # now works! log_scale = FALSE n=692; F=0.55028725698 after A3-A4 loop n= 3833 at convergence > ## Now that pbetaD94() can use 'log_scale', this converges at n=64974 (!): > pbetaD94(.99, a, b, ncp = 100, log_scale= TRUE, verbose=TRUE) -> pbL log_scale = TRUE n=692; F=0.55028725698 after A3-A4 loop n= 3833 at convergence > > all.equal(0.9999976, pb.) [1] TRUE > all.equal(0.9999976, pbL) [1] TRUE > > showProc.time() Time (user system elapsed): 0.07 0 0.07 > > ## Other ideas: > ## 1) in case Newton is not usable, be better than x' = (1+x)/2 {on right hand} : rather use Regula Falsi, or smart unirootR() ! > ## 2) in these cases, use rough estimates, e.g., a few steps of unirootR() ## with e.g., eps=1e-3 > > proc.time() user system elapsed 58.03 0.65 58.67