R Under development (unstable) (2024-07-29 r86934 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. > library(longmemo) > > capabilities() jpeg png tiff tcltk X11 aqua TRUE TRUE TRUE TRUE FALSE FALSE http/ftp sockets libxml fifo cledit iconv TRUE TRUE FALSE TRUE FALSE TRUE NLS Rprof profmem cairo ICU long.double TRUE TRUE TRUE TRUE TRUE TRUE libcurl TRUE > str(.Machine) List of 29 $ double.eps : num 2.22e-16 $ double.neg.eps : num 1.11e-16 $ double.xmin : num 2.23e-308 $ double.xmax : num 1.8e+308 $ double.base : int 2 $ double.digits : int 53 $ double.rounding : int 5 $ double.guard : int 0 $ double.ulp.digits : int -52 $ double.neg.ulp.digits : int -53 $ double.exponent : int 11 $ double.min.exp : int -1022 $ double.max.exp : int 1024 $ integer.max : int 2147483647 $ sizeof.long : int 4 $ sizeof.longlong : int 8 $ sizeof.longdouble : int 16 $ sizeof.pointer : int 8 $ sizeof.time_t : int 8 $ longdouble.eps : num 1.08e-19 $ longdouble.neg.eps : num 5.42e-20 $ longdouble.digits : int 64 $ longdouble.rounding : int 5 $ longdouble.guard : int 0 $ longdouble.ulp.digits : int -63 $ longdouble.neg.ulp.digits: int -64 $ longdouble.exponent : int 15 $ longdouble.min.exp : int -16382 $ longdouble.max.exp : int 16384 > str(.Platform) List of 8 $ OS.type : chr "windows" $ file.sep : chr "/" $ dynlib.ext: chr ".dll" $ GUI : chr "RTerm" $ endian : chr "little" $ pkgType : chr "win.binary" $ path.sep : chr ";" $ r_arch : chr "x64" > sessionInfo() # will change often.. but if we need the info, get it here R Under development (unstable) (2024-07-29 r86934 ucrt) Platform: x86_64-w64-mingw32/x64 Running under: Windows Server 2022 x64 (build 20348) Matrix products: default locale: [1] LC_COLLATE=C LC_CTYPE=German_Germany.utf8 [3] LC_MONETARY=C LC_NUMERIC=C [5] LC_TIME=C time zone: Europe/Berlin tzcode source: internal attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] longmemo_1.1-3 loaded via a namespace (and not attached): [1] compiler_4.5.0 > > ###-------------- Comparison of Paxson's approximation to B(): ----------------- > > B.o.specFGN <- function(lam, H, nsum=200) + B.specFGN(lam, H, k.approx=NA, nsum=nsum) > all.equal(B.specFGN (.ffreq(1000), 0.8), + B.o.specFGN(.ffreq(1000), 0.8)) [1] "Mean relative difference: 4.85353e-05" > ## 4.853 e-5 > > relErr <- function(target, current) { + xn <- mean(abs(target)) + stopifnot(is.finite(xn)) + if(xn == 0) { warning("target=0 -> using absolute average error"); xn <- 1 } + mean(abs(target-current)) / xn + } > > ## the relative error of the approxmation depends on H, not on n > ## and seems worst for H = 0.5 > H.s <- seq(0.5, 1, by= 1/32) > n.s <- c(25, 50, 100, 200, 500, 999, 2001) > p0 <- function(...) paste(..., sep="") > n.H <- function(H) p0("H=",formatC(H,digits=2,width=1)) > n.n <- function(n) p0("n=",as.character(n)) > rE <- matrix(NA, length(H.s), length(n.s), + dimnames=list(n.H(H.s), n.n(n.s))) > > for(H in H.s) { + ## the relative error of the approxmation depends on H, not on n + ## and seems worst for H = 0.5 + cat("H =", H,": ") + for(n in n.s) { + cat(".") + fr <- .ffreq(n) + rE[n.H(H), n.n(n)] <- + relErr(target = B.o.specFGN(fr, H, nsum=1000), + current= B.specFGN (fr, H)) + ## stopifnot(all.equal(B.specFGN (fr, H), + ## B.o.specFGN(fr, H, nsum=500), + ## tol = .25/n)) + } + cat("\n") + } H = 0.5 : ....... H = 0.53125 : ....... H = 0.5625 : ....... H = 0.59375 : ....... H = 0.625 : ....... H = 0.65625 : ....... H = 0.6875 : ....... H = 0.71875 : ....... H = 0.75 : ....... H = 0.78125 : ....... H = 0.8125 : ....... H = 0.84375 : ....... H = 0.875 : ....... H = 0.90625 : ....... H = 0.9375 : ....... H = 0.96875 : ....... H = 1 : ....... > > round(log10(rE), 2) n=25 n=50 n=100 n=200 n=500 n=999 n=2001 H=0.5 -3.37 -3.36 -3.36 -3.36 -3.36 -3.36 -3.36 H=0.53 -3.58 -3.57 -3.57 -3.57 -3.57 -3.57 -3.57 H=0.56 -3.79 -3.79 -3.78 -3.78 -3.78 -3.78 -3.78 H=0.59 -4.00 -3.99 -3.99 -3.99 -3.99 -3.99 -3.99 H=0.62 -4.20 -4.20 -4.20 -4.20 -4.19 -4.19 -4.19 H=0.66 -4.42 -4.41 -4.41 -4.41 -4.40 -4.40 -4.40 H=0.69 -4.66 -4.65 -4.64 -4.64 -4.64 -4.64 -4.64 H=0.72 -4.88 -4.86 -4.86 -4.86 -4.86 -4.86 -4.86 H=0.75 -4.99 -4.98 -4.97 -4.97 -4.96 -4.96 -4.96 H=0.78 -4.93 -4.94 -4.92 -4.91 -4.91 -4.91 -4.91 H=0.81 -4.72 -4.73 -4.72 -4.72 -4.72 -4.71 -4.71 H=0.84 -4.55 -4.57 -4.56 -4.56 -4.56 -4.56 -4.56 H=0.88 -4.43 -4.45 -4.44 -4.44 -4.44 -4.44 -4.44 H=0.91 -4.34 -4.35 -4.35 -4.35 -4.34 -4.34 -4.34 H=0.94 -4.26 -4.27 -4.27 -4.27 -4.27 -4.26 -4.26 H=0.97 -4.19 -4.21 -4.20 -4.20 -4.20 -4.20 -4.20 H=1 -4.14 -4.15 -4.15 -4.14 -4.14 -4.14 -4.14 > > cat('Time elapsed: ', proc.time(),'\n') # "statistics" Time elapsed: 12.7 0.59 13.28 NA NA > > proc.time() user system elapsed 12.70 0.59 13.28