R Under development (unstable) (2023-11-02 r85465 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > ### Testing numerical stability of h() {incl. Taylor around 0} > ### --------------------------- === > > library(DPQ) > ## --> pchisq.W(), h0(), h1(), h2() === h(), etc > > mult.fig <- sfsmisc::mult.fig > ## and (used via :: below) rrange <- sfsmisc::rrange > > ### TODO: Use many stopifnot() below <<<<<<<<<<<<< TODO >>> > > ## h() only defined in [0, 1] > > if(!dev.interactive(orNone=TRUE)) pdf("wienergerm-h-gnt.pdf") > > x <- seq(0, 1, length=101) > all.equal(h0(x), h1(x), tol= 5e-14)## TRUE (386 Linux) [1] TRUE > > ## not only in [0,1] : > plot(h, -4, 1, n= 2001) > abline(0, 1/6,col=2,lty=2); abline(h=c(0,1/2),v=0:1,col='gray', lty=3) > plot(h, 0, 1.1, n= 2001)## only in [0,1] Warning message: In log1p(-y.) : NaNs produced > abline(0, 1/6,col=2,lty=3); abline(h=c(0,1/2),v=0:1,col='gray',lty=3) > ## quite linear already: > > ## FIXME use h0,h1,h2 > plot(h, 0, 0.1, n= 2001);abline(0, 1/6,col=2,lty=3) > ## very linear: > plot(h, 0, 0.01, n= 2001);abline(0, 1/6,col=2,lty=3) > ## already breakdown for "h0"; ok for "h1" > plot(h, 0, 1e-4, n= 2001);abline(0, 1/6,col=2,lty=3) > ## noise even for "h1" (perfekt for "h2"): > plot(h, 0, 1e-7, n= 2001);abline(0, 1/6,col=2,lty=3) > ## Difference: looks very quadratic > x <- seq(0, 1e-7, length=2001) > plot(x, h2(x) - x/6, type='l', col=2)# look at y-scale! > > > ### Look close to 0 : > y <- 2^seq(-50,-11, length=2001) > y <- 2^seq(-24,-11, length=2001) > op <- mult.fig(2)$old.par > for(log in c('', "x")) { + plot(y,hnt(y,0)/y, type='l', ylim = 1/6+c(-1,1)*1e-3, log=log) + lines(y,hnt(y, 2)/y, col=3) + if(log=='')title("hnt[d](y)/y & hnt(y,2)/y {only terms k=0,1,2}") + else title("at y <- 2^ seq(-80,-10, length=2001)") + ## lines(y,hnt(y,10), col=2) + }; par(op) > > y <- 2^seq(-19,-11, length=2001) > op <- mult.fig(2)$old.par > for(log in c('', "x")) { + plot (y,hnt(y,0)/y - 1/6, type='l', ylim = sfsmisc::rrange(hnt(y,0)/y - 1/6), + log=log, xlab=NA) + lines(y,hnt(y,4)/y - 1/6, col=3) + if(log=='')title("hnt[d](y)/y -1/6 & hnt(y,4)/y - 1/6", xlab = "y") + else title("at y <- 2^ seq(-19,-11, length=2001)", xlab = "y [log scale]") + ## lines(y,hnt(y,10), col=2) + }; par(op) > > > ### Wiener germ approx of pchisq(x, df, ncp) -- Testing the Fortran / C code > ### ------------------ ------------ > ### ==> ./wienergerm-pchisq-tst.R > ### ======================= > > nu.lam.expr <- function(df, ncp) + ## a 'title' string + substitute(list(nu == df, lambda == ncp), + list(df = df, ncp = ncp)) > > ### ---> from now on checking C version only ==> pchisqW() <<<<< > ### -------------- --------- > > p.W.pchisq <- function(df, ncp, n = 513, x = NULL, do.log = TRUE, + cex = 0.5, s.col = "purple") + { + ## Purpose: Plot "Wiener germ" approx. pnchisq() problem around s=1 <==> x = df+ncp + ## ---------------------------------------------------------------------- ========== + ## Author: Martin Maechler, Date: 26 Jan 2004, 17:32 + + x0 <- df + ncp + if(is.null(x)) + x <- seq(0.5*min(df,ncp), 1.5*x0, length = n)## << better ? + + if(do.log) { + op <- mult.fig(2, main="wienergerm : problem at x ~= ncp+df (df 'small')", + marP = c(0,0,0,1.5), quiet = TRUE)$old.par + on.exit(par(op)) + } + clr <- c("red", adjustcolor(c("forestgreen","midnightblue"), 3/4, 1/6)) + for(log in c('x', if(do.log)'xy')) { + plot (x, pchisq (x, df=df, ncp=ncp), col = clr[1], log = log, cex = cex) + abline(h = 0.5, col = "gray", lty=3) + lines(x, pchisqW(x, df=df, ncp=ncp, var = "f"), col = clr[2], lwd=3) + lines(x, pchisqW(x, df=df, ncp=ncp, var = "s"), col = clr[3], lwd=5) + points(x0, pchisqW(x0, df=df, ncp=ncp), col = 1, pch=8, cex = 1.5) + axis(3, at=x0, label=quote(x[0] == nu + lambda)) + if(log == 'x') { + title(nu.lam.expr(df,ncp)) + legend(x[1], par('usr')[4],c("true","'f'","'s'"), + col=clr, pch= c(1,NA,NA), lty=c(NA,1,1), lwd=c(NA,3,5), + xjust=0, yjust=1.4, horiz=TRUE) + op <- par(new=TRUE) + sx <- sW(x, df=df, ncp=ncp)$s + plot(x, sx, type = 'l', col= s.col, log=log, axes=FALSE, ylab='') + x.mx <- 10^par("usr")[2] + segments(x0, 1, x.mx, 1, col=s.col, lty=2) + axis(4, col.axis = s.col, col = s.col) + mtext("s(x)",side=4, col = s.col, adj = 1.02) + } else { + title(sprintf('log = "%s"', log)) + } + } + + }# p.W.pchisq() > > p.W.pchisq(2 , 1) nonc_chi(*, variant=2) -> s = 1 ==> setting variant := 1 nonc_chi(*, variant=2) -> s = 1 ==> setting variant := 1 nonc_chi(*, variant=2) -> s = 1 ==> setting variant := 1 nonc_chi(*, variant=2) -> s = 1 ==> setting variant := 1 Warning message: In mtext(main, side = 3, outer = TRUE, line = line.main, cex = cex.main, : "quiet" is not a graphical parameter > p.W.pchisq(1 , 2) nonc_chi(*, variant=2) -> s = 1 ==> setting variant := 1 nonc_chi(*, variant=2) -> s = 1 ==> setting variant := 1 nonc_chi(*, variant=2) -> s = 1 ==> setting variant := 1 nonc_chi(*, variant=2) -> s = 1 ==> setting variant := 1 Warning message: In mtext(main, side = 3, outer = TRUE, line = line.main, cex = cex.main, : "quiet" is not a graphical parameter > p.W.pchisq(1 , 1) nonc_chi(*, variant=2) -> s = 1 ==> setting variant := 1 nonc_chi(*, variant=2) -> s = 1 ==> setting variant := 1 Warning message: In mtext(main, side = 3, outer = TRUE, line = line.main, cex = cex.main, : "quiet" is not a graphical parameter > p.W.pchisq(1 , 0.5)# here 'f' is quite different from 's' nonc_chi(*, variant=2) -> s = 1 ==> setting variant := 1 nonc_chi(*, variant=2) -> s = 1 ==> setting variant := 1 nonc_chi(*, variant=2) -> s = 1 ==> setting variant := 1 nonc_chi(*, variant=2) -> s = 1 ==> setting variant := 1 Warning message: In mtext(main, side = 3, outer = TRUE, line = line.main, cex = cex.main, : "quiet" is not a graphical parameter > p.W.pchisq(1 , 0.1) Warning message: In mtext(main, side = 3, outer = TRUE, line = line.main, cex = cex.main, : "quiet" is not a graphical parameter > ## df < 1 : becomes very bad (for s <= 1); not really good even further up > p.W.pchisq(.5, 0.1) Warning message: In mtext(main, side = 3, outer = TRUE, line = line.main, cex = cex.main, : "quiet" is not a graphical parameter > p.W.pchisq(.1, 0.2, x=2^seq(-50,3, length=2001)) nonc_chi(*, variant=2) -> s = 1 ==> setting variant := 1 nonc_chi(*, variant=2) -> s = 1 ==> setting variant := 1 Warning message: In mtext(main, side = 3, outer = TRUE, line = line.main, cex = cex.main, : "quiet" is not a graphical parameter > ## for larger x : "first" is clearly better than "second"! > p.W.pchisq(.1, 0.2, x=2^seq(-4,20,length=3001)) nonc_chi(*, variant=2) -> s = 1 ==> setting variant := 1 nonc_chi(*, variant=2) -> s = 1 ==> setting variant := 1 Warning message: In mtext(main, side = 3, outer = TRUE, line = line.main, cex = cex.main, : "quiet" is not a graphical parameter > > ## df larger --> (always ?) problem smaller > p.W.pchisq(4 , 1) nonc_chi(*, variant=2) -> s = 1 ==> setting variant := 1 nonc_chi(*, variant=2) -> s = 1 ==> setting variant := 1 Warning message: In mtext(main, side = 3, outer = TRUE, line = line.main, cex = cex.main, : "quiet" is not a graphical parameter > p.W.pchisq(3 , 2) Warning message: In mtext(main, side = 3, outer = TRUE, line = line.main, cex = cex.main, : "quiet" is not a graphical parameter > > p.W.pchisq(4 , 5) nonc_chi(*, variant=2) -> s = 1 ==> setting variant := 1 nonc_chi(*, variant=2) -> s = 1 ==> setting variant := 1 Warning message: In mtext(main, side = 3, outer = TRUE, line = line.main, cex = cex.main, : "quiet" is not a graphical parameter > p.W.pchisq(7 , 3) Warning message: In mtext(main, side = 3, outer = TRUE, line = line.main, cex = cex.main, : "quiet" is not a graphical parameter > p.W.pchisq(7 , 3, x = seq(8,100, length=501)) Warning message: In mtext(main, side = 3, outer = TRUE, line = line.main, cex = cex.main, : "quiet" is not a graphical parameter > > ## ncp >= 80 --- the now (2010++) important region: > p.W.pchisq(7 , 80, x = seq(8,100, length=501)) nonc_chi(*, variant=2) -> s = 1 ==> setting variant := 1 nonc_chi(*, variant=2) -> s = 1 ==> setting variant := 1 Warning message: In mtext(main, side = 3, outer = TRUE, line = line.main, cex = cex.main, : "quiet" is not a graphical parameter > curve(pchisqW(x, 7, 80), 84,89, n=1024); abline(h = 0.5, lty=2, col="gray50") > > ## how can I find the region where z < 0 ? > ## it's boundaries are where the final sign change should happen, not s <= 1! > > p.z.s <- function(df, ncp, n = 513, x = NULL, p.pchi = FALSE) + { + ## Purpose: Plot "Wiener germ" approx. pnchisq() problem + ## ---------------------------------------------------------------------- + ## Author: Martin Mächler, Date: 26 Jan 2004, 21:44 + + x0 <- df + ncp + if(is.null(x)) + x <- seq(0.6*mean(df,ncp), 1.25*x0, length = n) + else if(is.unsorted(x)) x <- sort(x) + + zfs <- cbind(z.f(x,df,ncp), + z.s(x,df,ncp)) + zmat <- cbind(z0 (x,df,ncp), zfs) + + ## Select interesting y-range + yl <- range(zmat) + ym <- min(zfs) + yl[1] <- pmax(5*ym, -0.2, yl[1]) + yl[2] <- pmin(0.5, 10*(-ym), yl[2]) + + ## select x-range from y-range + xl <- range(x[y.in <- apply(yl[1] <= zmat & zmat <= yl[2], 1, any)]) + ix <- xl[1] <= x & x <= xl[2] + ## and sub-set + x <- x[ix] + zmat <- zmat[ix, ,drop=FALSE] + zfs <- zfs [ix, ,drop=FALSE] + s <- sW(x,df,ncp)$s + + ## really interested in negative z's : + ineg <- range(which(apply(zfs < 0, 1, any)))# << first and last x-index + + if(p.pchi) { + op <- mult.fig(2, marP = c(0,0,0,1.5), quiet = TRUE)$old.par + on.exit(par(op)) + p.W.pchisq(df, ncp, do.log=FALSE) + } + matplot(s, zmat, type='l', ylim = yl, xlab = 's(.)', ylab = 'z(.)', + col = 2:4, lwd = c(1,2,1), lty = 1, main = nu.lam.expr(df,ncp)) + abline(h=0,lty=3) + + usr <- par('usr') + legend(s[1], usr[4],c("z0","z'f'","z's'"), + col=2:4, lwd= c(1,2,1), xjust=0, yjust=1.4, horiz=TRUE) + s0 <- s[ineg] + x0 <- x[ineg] + lab <- paste(formatC(s0), paste("x=",formatC(x0)), sep="\n") + op <- par(mgp=3:1) + axis(1, at= s0, labels=lab, pos= 0, lwd=2, col="midnight blue") + par(op) + list(x.range = x0, s.range = s0) + } ## p.z.s() > > > p.z.s(3,1, p.p= TRUE)# now shows 'z.s(s ~= 1)' problem, z.f() is fine $x.range [1] 3.0375 4.0000 $s.range [1] 0.7994565 1.0000000 > p.z.s(1,3, p.p= TRUE) nonc_chi(*, variant=2) -> s = 1 ==> setting variant := 1 $x.range [1] 2.679688 3.324219 $s.range [1] 0.7930244 0.8990963 > p.z.s(1,1, p.p= TRUE) nonc_chi(*, variant=2) -> s = 1 ==> setting variant := 1 $x.range [1] 0.7632812 1.4349609 $s.range [1] 0.5066187 0.7980605 > p.z.s(6,1, p.p= TRUE) $x.range [1] 6.154883 6.366113 $s.range [1] 0.8929273 0.9199634 > ## all these have df+ncp = 16.1 : > p.z.s(16, 0.1, p.p= TRUE) $x.range [1] 15.39697 15.47920 $s.range [1] 0.9565916 0.9616699 > p.z.s(14, 2.1, p.p= TRUE) $x.range [1] 15.29302 15.43042 $s.range [1] 0.9554311 0.9630524 > p.z.s(10, 6.1, p.p= TRUE) $x.range [1] 15.15918 15.35229 $s.range [1] 0.9571154 0.9660020 > p.z.s( 6,10.1, p.p= TRUE) $x.range [1] 15.02549 15.28369 $s.range [1] 0.9583184 0.9684597 > > ## (6, 10) for z.s(), (s ~= 1) problem: > curve(z.s(x,6,10), 16-.05, 16+ .05, ylim = c(0.014,0.018),n=2000) > ## "first" doesn't have the problem: > curve(z.f(x,6,10), add = TRUE, n = 2000, col=2) > > z.s(16.001, 6,10, verbose=TRUE) c(f= 4.33346, s= 1.00004, h= -6.41003e-06, qs=4.33331, z0=3.55211e-05, eta=-3.55025e-05, g=-5538.39, t1=16615.2, t2=-16615.3) [1] 0.01601603 > ## c(f= 4.33346, s= 1.00004, h= -6.41003e-06, qs=4.33331, z0=3.55211e-05, > ## eta=-3.55025e-05, g=-5538.39, t1=16615.2, t2=-16615.3) > ##[1] 0.01601603 > > p.zs1 <- function(df, ncp, + r.ex = c(-8, -5),# depends on (df, ncp) [not much though] + n = 513, signs = c(-1,1), type = 'l', + xlab = '', ylab = paste("z.s(x0 +/- exp(x)"), + main = nu.lam.expr(df,ncp) + ) + { + ## Purpose: problematic of z.s() when s --> 1 : + ## ---------------------------------------------------------------------- + ## Arguments: + ## ---------------------------------------------------------------------- + ## Author: Martin Maechler, Date: 27 Jan 2004, 13:33 + + x0 <- df + ncp + cx0 <- formatC(x0) + + xO <- pretty(exp(r.ex)) + xO <- xO[xO > 0] + xO <- c(xO[1]/2, xO) + ## x-limits : use "correct" on "+", artificially shifted left by Dx on "-" + dx <- diff(xl <- rx <- range(r.ex)) + e <- dx/16 + ## final size: dx + e + dx ; e ~= dx/10 + xl[1] <- xl[1] - (e+dx) + Dx <- dx+e ## = correction for left part + + at <- log(xO) + at <- list(2*rx[1]-e - at, at)# for (-1 , 1) + + xo <- seq(r.ex[1], r.ex[2], length = n) + + exo <- exp(xo) + + ## Compute everything -> y-layout + z <- as.list(signs) + for(is in seq(along=signs)) { + sig <- signs[is] + z[[is]] <- z.s(x0 + sig*exo, df, ncp) + } + yl <- sfsmisc::rrange(unlist(z), r = 2) # y-limits + plot(NA,NA, xlim = xl, ylim= yl, type = "n", xaxt = "n", + xlab=xlab, ylab=ylab, main=main) + uy <- par('usr')[3:4] + polygon(x= rx[1]+c(-e,-e,0,0), + y= uy[c(1:2,2:1)], col = 'gray70', xpd = FALSE) + mtext(paste(cx0, "+"), side = 1, cex = 1.5, adj = 0.5, line = 0.1) + op <- par(mgp=c(3,0.5,0)) + axis(3, at = pretty(rx), cex.axis = 0.8) + par(op) + mtext(paste("log( x -", cx0,")"), line= 1.2, adj = 0.75) + for(is in seq(along=signs)) { + sig <- signs[is] + lines(if(sig > 0) xo else rev(xo) - Dx, + y = z[[is]], type = type) + axis(1, at=at[[is]], labels = formatC(sig*xO)) + } + } > > p.zs1(3,1) > p.zs1(6,2) > p.zs1(6,10) > p.zs1(10,10) > p.zs1(10,12) > p.zs1(10,20)## oops! doesn't happen at "30" but slightly afterwards! > > p.zs1(10,200, c(-10,-3))# (two regions)! > > > ###-- Find improved z.s() --- need g(1-s) accurately ---------- > > tit.zs <- "pnchisq() Wiener germ z.'s' around s ~=1" > > curve(gnt(x, 0),-2, 1, n =1001, xlab = "u = 1 - s", main=tit.zs) > curve(gnt(x,10),-2, 1, n =1001, add=TRUE,col=2) > curve(gnt(x, 4),-2, 1, n =1001, add=TRUE,col=3) > curve(gnt(x, 2),-2, 1, n =1001, add=TRUE,col=4) > lines(0, gnt(0,2), type = "h", col="blue3", lty=3) > legend(-1,1.8, paste("g(u)", c("direct",paste(format(c(10,4,2)),"-term",sep=''))), + lty=1, col = 1:4) > > ### Look close to 0 : > y <- 2^seq(-50,-10, length=2001) > for(log in c('', "x")) { + plot(y,gnt(y,0), type='l', ylim = c(.499,.501), log=log) + lines(y,gnt(y, 2), col=3) + if(log=='')title("gnt[d](y) & gnt(y,2) {only terms k=0,1,2}") + else title("at y <- 2^ seq(-50,-10, length=2001)") + ## lines(y,gnt(y,10), col=2) + }; par(op) > > ## -- now look at relative *error* > > p.g <- function(to = 0.1, from = -to, Ns = c(1:4,10), cutoff = 0.2, + abs.do = FALSE, do.log = abs.do, leg.do = TRUE, + n.out = 513, main = tit.zs) + { + ## Purpose: + ## ---------------------------------------------------------------------- + ## Arguments: + ## ---------------------------------------------------------------------- + ## Author: Martin Maechler, Date: 29 Jan 2004, 13:51 + x <- seq(from,to, length=n.out) + g0 <- ifelse(abs(x) < cutoff, gnt(x,20), gnt(x,0)) + nN <- length(Ns) + gmat <- matrix(NA, n.out, nN) + for(j in seq(along=Ns)) { N <- Ns[j]; gmat[,j] <- gnt(x, N) } + y <- g0 / gmat - 1 + if(abs.do || do.log) + y <- abs(y) + if(do.log) { i0 <- y < 2*.Machine$double.eps; y[i0] <- NA } + matplot(x, y, type = 'l', ylim = sfsmisc::rrange(y), log = if(do.log)"y" else '', + xlab = '1-s', main = main, + ylab = paste("rel.Error g*(u) / g(u, n=",deparse(Ns),") - 1",sep='')) + yU <- par('usr')[3:4] ; if(do.log) yU <- 10^yU + if(leg.do) + legend("topleft",#0, if(abs.do) yU[2] else .7*yU[1]+.3*yU[2], + paste(format(Ns),"-term",sep=''), bty="n", + col = 1:nN, lty=1:nN, ncol = ceiling(nN/3)) + } > > p.g(1) > p.g(.1) > p.g(.001) > > > { op <- mult.fig(2)$old.par + p.g(.1) + p.g(.1, abs = TRUE, leg.do = FALSE) + par(op) } > > ## How bad is "direct", close to 0 : (looks worse than it is ?) > ## ## rrange > p.g(.001, Ns = 0)# +- 3e-5 > p.g(.004, Ns = 0)# +- 4e-7 > p.g(.020, Ns = 0)# +- 4e-9 > p.g(.040, Ns = 0)# +- 4e-10 > p.g(.1, Ns = 0)# +- 2e-11 > p.g(.2, Ns = 0)# +- 3e-12 > > > ### > epsC <- .Machine$double.eps > > ###===== Bound for simple two term (k = 0,1: linear) formula: =============== > (5/2*epsC) ^(1/2) ## 2.356 e-8 [1] 2.35608e-08 > (5/2*epsC/2)^(1/2) ## 1.666 e-8 [1] 1.666e-08 > ## slightly inside: > y <- 2.2 * 10^seq(-20,-8, length=1001) > > ## Absolute error: > summary(ae <- gnt(y,10) - gnt(y,1)) Min. 1st Qu. Median Mean 3rd Qu. Max. -1.110e-16 0.000e+00 0.000e+00 -2.118e-17 0.000e+00 1.110e-16 > ##- Min. 1st Qu. Median Mean 3rd Qu. Max. > ##- -1.110e-16 0.000e+00 0.000e+00 -2.118e-17 0.000e+00 1.110e-16 > > ## relative error: > summary(re <- 1 - gnt(y,1) / gnt(y,10)) Min. 1st Qu. Median Mean 3rd Qu. Max. -2.220e-16 0.000e+00 0.000e+00 -4.237e-17 0.000e+00 2.220e-16 > ##- Min. 1st Qu. Median Mean 3rd Qu. Max. > ##- -2.220e-16 0.000e+00 0.000e+00 -4.237e-17 0.000e+00 2.220e-16 > are <- abs(re) > > plot(y,are, log = 'x') > plot(y,are, log = 'x', xlim=c(1e-10,max(y))) > plot(y,are, log = 'xy', xlim=c(1e-10,max(y))) Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 782 y values <= 0 omitted from logarithmic plot > plot(y,are, log = '', xlim=c(1e-10,max(y))) > ##--> I should rather use epsC / 2 (makes sense from what we know!) > > ###===== Bound for three term (k = 0,1,2) formula: ==================== > (7/2*epsC) ^(1/3)## 9.19 e-6 [1] 9.193963e-06 > (7/2*epsC/2)^(1/3)## 7.30 e-6 [1] 7.297253e-06 > > ## slightly inside: > y <- 9* 10^seq(-20,-6, length=1001) > > ## Absolute error: > range(ae <- gnt(y,10) - gnt(y,2))## 0 1.11e-16 [1] 0.000000e+00 1.110223e-16 > > ## relative error: > range(re <- 1 - gnt(y,2) / gnt(y,10))## 0 2.22e-16 [1] 0.000000e+00 2.220446e-16 > > plot(y,re, log = 'x') > plot(y,re, log = 'x', xlim=c(1e-6,max(y))) > plot(y,re, log = 'xy', xlim=c(1e-6,max(y))) Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 986 y values <= 0 omitted from logarithmic plot > plot(y,re, log = '', xlim=c(1e-6,max(y))) > ###--> I should rather use epsC / 2 (makes sense from what we know!) > ### or even a bit less ..: |y| < 1e-5 > > > > ###===== Bound for four term (k = 0,1,2,3) formula: ==================== > (14/3*epsC) ^(1/4) ## 0.000179 [1] 0.0001794162 > (14/3*epsC/2)^(1/4) ## 0.000151 [1] 0.0001508705 > > ## slightly _out_side: > y <- 1.8* 10^seq(-10,-4, length=1001) > ## Absolute error: > range(ae <- gnt(y,10) - gnt(y,3)) # 0. 1.11e-16 [1] 0.000000e+00 1.110223e-16 > ## relative error (very similar) > range(re <- 1 - gnt(y,3) / gnt(y,10))# 0. 2.22e-16 [1] 0.000000e+00 2.220446e-16 > > > ###===== Bound for five term (k = 0,1,2,3,4) formula: ==================== > (6*epsC )^(1/5) ## 0.00106 [1] 0.001059054 > (6*epsC/2)^(1/5) ## 0.000922 [1] 0.0009219605 > > ## slightly inside: > y <- 1* 10^seq(-10,-3, length=1001) > ## Absolute error: > range(ae <- gnt(y,10) - gnt(y,4)) # 0 1.11e-16 [1] 0.000000e+00 1.110223e-16 > ## relative error (very similar) > range(re <- 1 - gnt(y,4) / gnt(y,10))# 0 2.22e-16 [1] 0.000000e+00 2.220446e-16 > > > ###===== Bound for 6 term (k = 0,1,2,3,4,5) formula: ==================== > (15/2*epsC )^(1/6) ## 0.00344 [1] 0.003442841 > > ## slightly inside: > y <- 3.4* 10^seq(-10,-3, length=1001) > ## Absolute error: > range(ae <- gnt(y,10) - gnt(y,5)) # 0 1.11e-16 [1] 0.000000e+00 1.110223e-16 > ## relative error (very similar) > range(re <- 1 - gnt(y,5) / gnt(y,10))# 0 2.22e-16 [1] 0.000000e+00 2.220446e-16 > > > ###===== Bound for seven term (k = 0,1,2,3,4,5,6) formula: ==================== > (55/6*epsC )^(1/7) ## 0.00797 [1] 0.00796592 > ## slightly outside: > y <- 8* 10^seq(-10,-3, length=1001) > ## Absolute error: > range(ae <- gnt(y,10) - gnt(y,6)) # 0 2.22-16 [1] 0.000000e+00 2.220446e-16 > ## relative error (very similar) > range(re <- 1 - gnt(y,6) / gnt(y,10))# 0 4.44e-16 [1] 0.000000e+00 4.440892e-16 > > ##==> MM{2014}: I think I settled for g2() in 2004, > ## -------- which goes up to 4-term ==== because it is good enough. > > > > ###===== former ../tests/wienergerm-accuracy.R ========================================== > ### Orig. = R script /u/maechler/R/MM/NUMERICS/dpq-functions/wienergerm-accuracy.R > > if(!dev.interactive(orNone=TRUE)) pdf("wienergerm-accuracy.pdf") > > nuS <- c(1:3, 10, 20, 30,50, 100, 200, 1000, 5000, 1e4) > ncpS <- c(0, .1, 1, 3, 10, 30, 100, 300, 1000, 10000, 1e5) > facs <- c(outer(c(1,2,5), 10^c(-1:3))) > r <- array(NA, dim = c(length(facs), 3, length(nuS), length(ncpS)), + dimnames= list(x=NULL, c("p.","pW","pP"), + df = formatC(nuS), ncp= formatC(ncpS))) > for(df in nuS) { + for(ncp in ncpS) { + m <- df + ncp + s <- sqrt(2*(df + 2*ncp)) + x <- m + s * facs + r[,, formatC(df), formatC(ncp)] <- + cbind(pchisq (x, df=df, ncp=ncp), + pchisqW (x, df=df, ncp=ncp), + pnchisqPearson(x, df=df, ncp=ncp)) + } + cat("done df=",formatC(df),"\n") + } done df= 1 done df= 2 done df= 3 done df= 10 done df= 20 done df= 30 done df= 50 done df= 100 done df= 200 done df= 1000 done df= 5000 done df= 1e+04 > > p.pchisq <- function(x, df, ncp, log = "x", + relErr = FALSE, diff = FALSE) + { + ## Purpose: + ## ---------------------------------------------------------------------- + ## Arguments: + ## ---------------------------------------------------------------------- + ## Author: Martin Maechler, Date: 18 Mar 2004, 17:04 + + Px <- pchisq(x, df=df, ncp=ncp) + pW <- pchisqW(x, df=df, ncp=ncp) # "Wiener2" + pW1 <- pchisqW(x, df=df, ncp=ncp, variant = 'f') + pP1 <- pnchisqPearson(x, df=df, ncp=ncp) + pP2 <- pnchisqPatnaik(x, df=df, ncp=ncp) + Pmat <- cbind(pW,pW1,pP1,pP2) + + tit <- paste("pchisq[Appr](x, df=",formatC(df),", ncp=",formatC(ncp),")") + if(diff) + tit <- paste(tit," - pchisq[R](.)") + else if( relErr) + tit <- paste("rel.Error of", tit," to pchisq[R](.)") + log.y <- length(grep("y", log)) + col <- 2:5 + lty <- 1:4 + noR <- diff || relErr # subtract from R's pchisq() -- TODO: rather use Rmpfr-ified pnchisq()! + if(noR) { + Pmat <- Pmat - Px + if(relErr) + Pmat <- abs(Pmat / Px) + else if(log.y) { + Pmat <- abs(Pmat) + if(log.y) tit <- paste("|", tit, "|") + } + } + else { + Pmat <- cbind(Px, Pmat) + col <- c(1,col) + lty <- c(1,lty) + } + if(log.y) op <- options(warn=-1)# no "log warning" + matplot (x, Pmat, log=log, type='l', main = tit, + ylab = '', ## ylim = rrange(Pmat, range = 0.1), + lty = lty, col = col) + if(log.y) op <- options(op) + u <- par("usr") + log.x <- length(grep("x", log)) + if(log.x) u[1:2] <- 10^u[1:2] + if(log.y) u[3:4] <- 10^u[3:4] + legend(u[1],u[4], + c(if(!noR)"R", "Wiener2", "Wiener1", "Pearson","Patnaik"), + lty = lty, col = col) + } > > ### all these for (0.1, 0.1): > x <- 2^seq(-10,3, length=501) > p.pchisq(x, df=0.1, ncp=0.1) > p.pchisq(2^seq(-12,6, length=2001), df=0.1, ncp=0.1) > p.pchisq(2^seq(-1, 6, length=2001), df=0.1, ncp=0.1) > p.pchisq(2^seq(-1, 6, length=2001), df=0.1, ncp=0.1, diff = TRUE) > > p.pchisq(2^seq(-1, 10, length=2001), df=0.1, ncp=0.1, diff = TRUE) > p.pchisq(2^seq(-1, 10, length=2001), df=0.1, ncp=0.1, diff = TRUE,log="xy") > p.pchisq(2^seq(-5, 8, length=2001), df=0.1, ncp=0.1, diff = TRUE,log="xy") > p.pchisq(2^seq(-5, 8, length=2001), df=0.1, ncp=0.1, relE = TRUE,log="xy") > > ## (1, 1) > p.pchisq(2^seq(-14,8, length=2001), df=1, ncp= 1) > p.pchisq(2^seq(-14,8, length=2001), df=1, ncp= 1,diff=TRUE) > p.pchisq(2^seq(-20,9, length=2001), df=1, ncp= 1,diff=TRUE,log="xy") > p.pchisq(2^seq(-20,9, length=2001), df=1, ncp= 1,relE=TRUE,log="xy") > > ## (1, 10) > p.pchisq(2^seq(-14,8, length=2001), df=1, ncp= 10,diff=TRUE,log="xy") > ## (10, 1) > p.pchisq(2^seq(-14,8, length=2001), df=10, ncp= 1,diff=TRUE,log="xy") > > ## (10, 10) > p.pchisq(2^seq(-14,8, length=2001), df=10, ncp= 10,diff=TRUE,log="xy") > p.pchisq(2^seq(-1,10, length=2001), df=10, ncp= 10,diff=TRUE,log="xy") > > ## (10, 100) > p.pchisq(2^seq(-14,8, length=2001), df=10, ncp= 100,diff=TRUE,log="xy") > p.pchisq(2^seq(-1,10, length=2001), df=10, ncp= 100,diff=TRUE,log="xy") > ## "good" news: Pearson is "nice" when Wiener fails > > ## (10, 1e4) > p.pchisq(seq(9000, 11000, length=2001), df=10, ncp= 1e4, log='') nonc_chi(*, variant=2) -> s = 1 ==> setting variant := 1 > p.pchisq(seq(9000, 11000, length=2001), df=10, ncp= 1e4, diff=TRUE) nonc_chi(*, variant=2) -> s = 1 ==> setting variant := 1 > p.pchisq(seq(9000, 11000, length=2001), df=10, ncp= 1e4, diff=TRUE, log='y') nonc_chi(*, variant=2) -> s = 1 ==> setting variant := 1 > p.pchisq(seq(9000, 11000, length=2001), df=10, ncp= 1e4, relE=TRUE, log='y') nonc_chi(*, variant=2) -> s = 1 ==> setting variant := 1 > > ## (1000, 1) > p.pchisq(seq(800, 1200, length=2001), df=1000, ncp= 1, log='') > > p.pchisq(seq(800, 1200, length=2001), df=1000, ncp= 1, diff=TRUE, log='y') > p.pchisq(seq(800, 1200, length=2001), df=1000, ncp= 1, relE=TRUE, log='y') > ## interesting: Pearson really better than Wiener in a whole region > ## x in [900, 1100] ... > > proc.time() user system elapsed 7.67 0.29 7.92