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. > ### originally was ~/R/MM/NUMERICS/dpq-functions/pnt-precision-problem.R > ### - - - - - - - - - - - - - - - - - - - - - - - - - - - > ##==> see also ./t-nonc-tst.R > ## ============ > library(DPQ) > > stopifnot(exprs = { + require(graphics) + require(sfsmisc) # lseq(), p.m(), mult.fig() + }) Loading required package: sfsmisc > > source(system.file(package="DPQ", "test-tools.R", + mustWork=TRUE))# ../inst/test-tools.R > ## => showProc.time(), ... 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("pnt-precision-2.pdf") > op <- options(nwarnings = 1e5) > > x <- 10^seq(2,12, by= .25) > px <- pt(x, df= 0.9, ncp = .01, + lower.tail=FALSE, log=TRUE) ## 16 warnings > ## full precision was not achieved in 'pnt' > plot(x,px, log="x", type="o") #- show catastrophic behavior > ## R-devel 2008-11-04: slope is kept, but we have a jump still > ## > ## extend x --> "the flattening" (at log(P[ ]) ~= -32) happens still > curve(pt(x, df= 0.9, ncp = .01, lower.tail=FALSE, log=TRUE), + 100, 1e20, log="x") > > xs <- lseq(100, 1e6, len = 101) > pxs <- pt(xs, df= 0.9, ncp = .01, lower.tail=FALSE, log=TRUE) > coef(lm(pxs ~ log(xs))) (Intercept) log(xs) -1.1484344 -0.8999986 > ## (Intercept) log(x) > ## -1.1484345 -0.8999986 > ## But (central) pt() has now such problem: > x <- 10^seq(2,17, by=1/4) > mp <- cbind(x, + pnt = pt(x, df= 0.9, ncp = .01, lower.tail=FALSE, log=TRUE), + pt = pt(x, df= 0.9, lower.tail=FALSE, log=TRUE)) > mp x pnt pt [1,] 1.000000e+02 -5.293100 -5.305215 [2,] 1.778279e+02 -5.811164 -5.823278 [3,] 3.162278e+02 -6.329240 -6.341354 [4,] 5.623413e+02 -6.847319 -6.859434 [5,] 1.000000e+03 -7.365400 -7.377515 [6,] 1.778279e+03 -7.883482 -7.895596 [7,] 3.162278e+03 -8.401564 -8.413678 [8,] 5.623413e+03 -8.919645 -8.931760 [9,] 1.000000e+04 -9.437727 -9.449841 [10,] 1.778279e+04 -9.955808 -9.967923 [11,] 3.162278e+04 -10.473890 -10.486005 [12,] 5.623413e+04 -10.991972 -11.004086 [13,] 1.000000e+05 -11.510053 -11.522168 [14,] 1.778279e+05 -12.028134 -12.040250 [15,] 3.162278e+05 -12.546221 -12.558331 [16,] 5.623413e+05 -13.064280 -13.076413 [17,] 1.000000e+06 -13.582351 -13.594494 [18,] 1.778279e+06 -14.100547 -14.112576 [19,] 3.162278e+06 -14.618350 -14.630658 [20,] 5.623413e+06 -15.137231 -15.148739 [21,] 1.000000e+07 -15.649607 -15.666821 [22,] 1.778279e+07 -16.183938 -16.184903 [23,] 3.162278e+07 -16.696748 -16.702984 [24,] 5.623413e+07 -17.139016 -17.221066 [25,] 1.000000e+08 -22.146393 -17.739148 [26,] 1.778279e+08 -22.664445 -18.257229 [27,] 3.162278e+08 -23.182477 -18.775311 [28,] 5.623413e+08 -23.700477 -19.293393 [29,] 1.000000e+09 -24.218417 -19.811474 [30,] 1.778279e+09 -24.736265 -20.329556 [31,] 3.162278e+09 -25.253953 -20.847638 [32,] 5.623413e+09 -25.771382 -21.365719 [33,] 1.000000e+10 -26.288347 -21.883801 [34,] 1.778279e+10 -26.804580 -22.401882 [35,] 3.162278e+10 -27.319527 -22.919964 [36,] 5.623413e+10 -27.832442 -23.438046 [37,] 1.000000e+11 -28.341775 -23.956127 [38,] 1.778279e+11 -28.845470 -24.474209 [39,] 3.162278e+11 -29.339852 -24.992291 [40,] 5.623413e+11 -29.820086 -25.510372 [41,] 1.000000e+12 -30.276896 -26.028454 [42,] 1.778279e+12 -30.698930 -26.546536 [43,] 3.162278e+12 -31.073840 -27.064617 [44,] 5.623413e+12 -31.389693 -27.582699 [45,] 1.000000e+13 -31.643050 -28.100781 [46,] 1.778279e+13 -31.831526 -28.618862 [47,] 3.162278e+13 -31.957677 -29.136944 [48,] 5.623413e+13 -32.045453 -29.655025 [49,] 1.000000e+14 -32.102072 -30.173107 [50,] 1.778279e+14 -32.141681 -30.691189 [51,] 3.162278e+14 -32.162090 -31.209270 [52,] 5.623413e+14 -32.172452 -31.727352 [53,] 1.000000e+15 -32.182924 -32.245434 [54,] 1.778279e+15 -32.182924 -32.763515 [55,] 3.162278e+15 -32.182924 -33.281597 [56,] 5.623413e+15 -32.182924 -33.799679 [57,] 1.000000e+16 -32.193506 -34.317760 [58,] 1.778279e+16 -32.193506 -34.835842 [59,] 3.162278e+16 -32.193506 -35.353924 [60,] 5.623413e+16 -32.193506 -35.872005 [61,] 1.000000e+17 -32.193506 -36.390087 > p.m(mp, type="l", log="x", xlab="x", main = "pt(x, df = 0.9, lower = F, log=TRUE)") > legend("topright", c("ncp = 0.01", "central"), lty=1:2, col=1:2, bty="n") > > ## even closer: (ncp ~= 0) : --> even bigger problem : pnt(..) becomes -Inf, > ## (even in R-devel 2008-11-04) > m2 <- cbind(x, + pnt = pt(x, df=.9, ncp=1e-4, lower.tail=FALSE, log=TRUE), + pnt = pt(x, df=.9, ncp=1e-10, lower.tail=FALSE, log=TRUE), + pt = pt(x, df=.9, lower.tail=FALSE, log=TRUE)) > m2 x pnt pnt pt [1,] 1.000000e+02 -5.305093 -5.305215 -5.305215 [2,] 1.778279e+02 -5.823157 -5.823278 -5.823278 [3,] 3.162278e+02 -6.341233 -6.341354 -6.341354 [4,] 5.623413e+02 -6.859312 -6.859434 -6.859434 [5,] 1.000000e+03 -7.377394 -7.377515 -7.377515 [6,] 1.778279e+03 -7.895475 -7.895596 -7.895596 [7,] 3.162278e+03 -8.413557 -8.413678 -8.413678 [8,] 5.623413e+03 -8.931638 -8.931760 -8.931760 [9,] 1.000000e+04 -9.449720 -9.449841 -9.449841 [10,] 1.778279e+04 -9.967802 -9.967923 -9.967923 [11,] 3.162278e+04 -10.485883 -10.486005 -10.486005 [12,] 5.623413e+04 -11.003965 -11.004086 -11.004086 [13,] 1.000000e+05 -11.522046 -11.522168 -11.522168 [14,] 1.778279e+05 -12.040127 -12.040248 -12.040250 [15,] 3.162278e+05 -12.558214 -12.558336 -12.558331 [16,] 5.623413e+05 -13.076273 -13.076394 -13.076413 [17,] 1.000000e+06 -13.594344 -13.594466 -13.594494 [18,] 1.778279e+06 -14.112541 -14.112663 -14.112576 [19,] 3.162278e+06 -14.630341 -14.630462 -14.630658 [20,] 5.623413e+06 -15.149232 -15.149353 -15.148739 [21,] 1.000000e+07 -15.661539 -15.661659 -15.666821 [22,] 1.778279e+07 -16.196066 -16.196189 -16.184903 [23,] 3.162278e+07 -16.708812 -16.708934 -16.702984 [24,] 5.623413e+07 -17.150194 -17.150307 -17.221066 [25,] 1.000000e+08 -26.755241 -Inf -17.739148 [26,] 1.778279e+08 -27.273370 -Inf -18.257229 [27,] 3.162278e+08 -27.791338 -Inf -18.775311 [28,] 5.623413e+08 -28.309532 -Inf -19.293393 [29,] 1.000000e+09 -28.827678 -Inf -19.811474 [30,] 1.778279e+09 -29.345385 -Inf -20.329556 [31,] 3.162278e+09 -29.863637 -Inf -20.847638 [32,] 5.623413e+09 -30.380693 -Inf -21.365719 [33,] 1.000000e+10 -30.899070 -Inf -21.883801 [34,] 1.778279e+10 -31.418681 -Inf -22.401882 [35,] 3.162278e+10 -31.932780 -Inf -22.919964 [36,] 5.623413e+10 -32.446341 -Inf -23.438046 [37,] 1.000000e+11 -32.975600 -Inf -23.956127 [38,] 1.778279e+11 -33.478704 -Inf -24.474209 [39,] 3.162278e+11 -33.964212 -Inf -24.992291 [40,] 5.623413e+11 -34.539576 -Inf -25.510372 [41,] 1.000000e+12 -34.945041 -Inf -26.028454 [42,] 1.778279e+12 -35.638188 -Inf -26.546536 [43,] 3.162278e+12 -36.043653 -Inf -27.064617 [44,] 5.623413e+12 -36.736801 -Inf -27.582699 [45,] 1.000000e+13 -36.736801 -Inf -28.100781 [46,] 1.778279e+13 -36.736801 -Inf -28.618862 [47,] 3.162278e+13 -Inf -Inf -29.136944 [48,] 5.623413e+13 -Inf -Inf -29.655025 [49,] 1.000000e+14 -Inf -Inf -30.173107 [50,] 1.778279e+14 -Inf -Inf -30.691189 [51,] 3.162278e+14 -Inf -Inf -31.209270 [52,] 5.623413e+14 -Inf -Inf -31.727352 [53,] 1.000000e+15 -Inf -Inf -32.245434 [54,] 1.778279e+15 -Inf -Inf -32.763515 [55,] 3.162278e+15 -Inf -Inf -33.281597 [56,] 5.623413e+15 -Inf -Inf -33.799679 [57,] 1.000000e+16 -Inf -Inf -34.317760 [58,] 1.778279e+16 -Inf -Inf -34.835842 [59,] 3.162278e+16 -Inf -Inf -35.353924 [60,] 5.623413e+16 -Inf -Inf -35.872005 [61,] 1.000000e+17 -Inf -Inf -36.390087 > p.m(m2, type="l", log="x", xlab="x", main = "pt(x, df = 0.9, lower = F, log=TRUE)") > legend("topright", c("ncp = 1e-4", "ncp = 1e-10", "central"), lty=1:3, col=1:3, bty="n") > > ### Note: The non-central F ( <==> non-central beta ) is sometimes *WORSE*: > p.t.vs.F <- function(df, ncp, x1 = 1, x2= 1e5, nout = 201, col = c("blue","red")) + { + ## Purpose: Tail - Comparison of non-central t & non-central F + ## Abramowitz & Stegun 26.6.19, p.947 + + ###>> This "non-sense" non-central F <==> "two-tail" non-central t + ###>> i.e., there's NO direct representation *unless* ncp = 0 + + ## ---------------------------------------------------------------------- + ## Arguments: + ## ---------------------------------------------------------------------- + ## Author: Martin Maechler, Date: 1 Nov 2008, 17:12 + x <- lseq(x1, x2, length=nout) + m <- cbind(pt=pt(x, df=df, ncp= ncp , log=TRUE,lower=FALSE), + pf=pf(x^2, df1=1, df2=df, ncp= ncp^2, log=TRUE,lower=FALSE)) + matplot(x, m, type="l", log="x", ylab = "log(1 - F*(x))", lwd = 2, col=col, + main = "upper tail prob. p*(......, log=TRUE, lower.tail=FALSE)") + legend("topright", + c(sprintf("pt(x, df = %g, ncp=%g )", df,ncp), + sprintf("pf(x^2, df1=1, df2= %g, ncp=%g ^2)",df,ncp)), + lty=1:2, col=col, lwd = 2, inset=.01) + mtext(R.version.string, side=1, line=-1, adj=0.01, cex = 0.6) + invisible(cbind(x=x, m)) + } > > p.t.vs.F(3.14, 2) > p.t.vs.F( 2, 20,, 1e9) > p.t.vs.F(0.2, 20,, 1e50)# pnt()_R-devel has small jump; 2.8.0 big jumps to horizontal > > ## but here (df << ncp), pf() {ie pnbeta} is clearly better}: > p.t.vs.F(2, 200, 10,1e6) > p.t.vs.F(2, 200, 10,1e9)# breaks down eventually as well > > ## Small ncp, *however* one of the two is WAY wrong ! -- "SYSTEMATIC ERROR" : > p.t.vs.F(3, 0.1) ## {and pf() breaks earlier for large x} There were 61 warnings (use warnings() to see them) > p.t.vs.F(3, 0.5) ## still {slightly} There were 50 warnings (use warnings() to see them) > p.t.vs.F(3, 1) ## no visible problem > ## but here, for x -> 0: pf() is systematically larger : > p.t.vs.F(3, 1, 1e-2,100) > > showProc.time() Time (user system elapsed): 0.19 0.02 0.2 > > ### --------------------------------------------------------------------------------- > > ## -> pntR() is R-level pnt() function > pt (1e10, df=.9, ncp=1e-10, lower.tail=FALSE, log=TRUE) # >> -Inf [1] -Inf > pntR(1e10, df=.9, ncp=1e-10, lower.tail=FALSE, log=TRUE) # " pnt(t= 1e+10 , df= 0.9 , delta= 1e-10 ) ==> x= 1 : p= 0.5 1 iter.: tnc{sum} = 0.50000000004, tnc+pnorm(-del) = 1, lower.t = 0 [1] -Inf > ## "finis: precision NOT reached" > > ## back to less extreme ncp: > x <- 2^(0:40) > px <- pt(x, df = 1., ncp = 1, lower.tail = FALSE, log = TRUE)## 14 warnings > plot(x,px, log="x", type="o") #- show "jump" {just small kink: R-devel 25.10.09} > > ## also the left tail (!) --- this is now (R-devel 25.Oct.09) worse > px <- pt(-x, df = 1.2, ncp = 1, log = TRUE) > plot(x, px, log="x", type="o") #- jump and wrong > > ## larger df: no jump, still wrong > px <- pt(x, df = 2, ncp = 1, lower.tail = FALSE, log = TRUE) ## 23 warnings > plot(x, px, log="x", type="o") #- show bad behavior > ## does my one approximation help here ? > px. <- pntJW39(x, df = 2, ncp = 1, lower.tail = FALSE, log = TRUE) > lines(x,px., col = 2, lwd=2) > ## no! It's worse for large x! > > px <- pt(x, df = 2, ncp = 2, lower.tail = FALSE, log = TRUE) ## 23 warnings > plot(x,px, log="x", type="o") #- show bad behavior > > ## The Mathematica example > ## LogLogPlot[ 1 - CDF[NoncentralStudentTDistribution[3, 5], x], {x, 1, 10^6}] > ## takes many minutes of computing !!!! > px <- pt(x, df = 3, ncp = 5, lower.tail = FALSE, log = TRUE) > plot(x, px, log="x", type="o") #- show bad behavior > ## equivalent {different x-labels}: > plot(log(x), px, type="o") > > ## BTW: Very slow convergence here {not in outer tail} -- can clearly be improved > pntR(20, df=0.9, ncp=30, lower.tail=FALSE, log=TRUE, errmax=1e-15) pnt(t= 20 , df= 0.9 , delta= 30 ) ==> x= 0.9977551 : p= 1.846942e-196 621 iter.: tnc{sum} = 0.136209733068, tnc+pnorm(-del) = 0.136209733068, lower.t = 0 [1] -0.1464253 > pntR(20, df= 2, ncp=30, lower.tail=FALSE, log=TRUE, errmax=1e-15) pnt(t= 20 , df= 2 , delta= 30 ) ==> x= 0.9950249 : p= 1.846942e-196 619 iter.: tnc{sum} = 0.106320225789, tnc+pnorm(-del) = 0.106320225789, lower.t = 0 [1] -0.1124078 > > ### "Solve this problem": Find tail formula empirically > ### ---------------------------------------- > ### log(1 - P(x)) ~= alpha - beta * log(x) > ### ---------------------------------------- > > ## Empirically: alpha = g(ncp) ; beta = h(df) << see below: beta == df ( = nu) !! > > summary(lm(px ~ log(x), subset=x > 7 & x < 35000))# seems clear, R^2 = 0.9999 Call: lm(formula = px ~ log(x), subset = x > 7 & x < 35000) Residuals: Min 1Q Median 3Q Max -0.25113 -0.03421 0.02592 0.05383 0.09301 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 5.03342 0.06876 73.2 3.84e-16 *** log(x) -2.96800 0.01018 -291.6 < 2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.09517 on 11 degrees of freedom Multiple R-squared: 0.9999, Adjusted R-squared: 0.9999 F-statistic: 8.504e+04 on 1 and 11 DF, p-value: < 2.2e-16 > > showProc.time() Time (user system elapsed): 0.15 0 0.16 > > ptRTailAsymp <- function(df,ncp, f.x1 = 1e5, f.x0 = 1, nx = 1000, + do.plot=FALSE, verbose = do.plot, + ask = do.plot) + { + ## Purpose: Right tail asymptotic of pt_noncentral() + ## ---------------------------------------------------------------------- + ## Arguments: + ## ---------------------------------------------------------------------- + ## Author: Martin Maechler, Date: 23 Oct 2008, 15:15 + stopifnot(length(df) == 1, is.numeric(df), df > 0) + stopifnot(length(ncp) == 1, is.numeric(ncp), ncp >= 0) + + n.x0 <- round(.25 * nx) # at this index, we anchor the robust line + ## Determine an x-range, where pt(*, lower.tail=FALSE, log=TRUE) + ## is nearly linear + alp <- .999 + x0. <- qt(alp, df=df, ncp=ncp)##<< it's absurd to use this here, + ## -- but unfortunately, qtAppr() is too unreliable here + it <- 1 + while(!is.finite(x0.) && it < 20) { + it <- it+1 + alp <- 1 - 2*(1-alp) + x0. <- qt(alp, df=df, ncp=ncp) + } + stopifnot(is.finite(x0.)) + x0 <- x0.* f.x0 + x1 <- x0 * f.x1 + it <- 1; converged <- FALSE + if(do.plot && ask) { op <- par(ask=TRUE); on.exit(par(op)) } + repeat { + xx <- lseq(x0, x1, len = nx) + lx <- log(xx) + px <- pt(xx, df=df, ncp=ncp, lower.tail=FALSE, log=TRUE) + rob.slope <- median(diff(px), na.rm=TRUE) /mean(diff(lx)) + Out <- !is.finite(px) + SSY <- sum((px[!Out] - mean(px[!Out]))^2) + if(do.plot) { + r <- lm.fit(cbind(1, lx[!Out]), px[!Out]) + if(it == 1) { ## first plot + plot(px ~ lx, xlab = quote(log(x)), + xlim = extendrange(r=log(c(x0,x1)), f=0.2), + ylim = extendrange(px, f = 0.25), + main = sprintf("pt(exp(.), df=%g, ncp=%g, lower.tail=FALSE, log.p=TRUE)", + df, ncp)) + mtext(sprintf("Search for [x0, x1] s.t. pt(*, log=TRUE) is *linear* in it; starting @ [%g, %g]", + log(x0), log(x1))) + } else { # it >= 2 ... subsequent iterations + lines(px ~ lx) + lines(lx[!Out], r$fitted, col= "green2") + abline(a = px[n.x0] - rob.slope*lx[n.x0], + b = rob.slope, col = "tomato", lwd=2) + } + } + rob.fitted <- px[n.x0] + rob.slope*(lx - lx[n.x0]) + rob.resid <- px - rob.fitted + + R2.rob <- 1 - sum(rob.resid^2) / SSY + cat(sprintf("robust R^2 = %12.8f\n", R2.rob)) + if(R2.rob < 0.99) { + x1 <- x1/10 + if(chg0 <- x0 > x1/2) x0 <- x0/2 + message(c(sprintf("'robust' fit not ok\n ==> x1 := x1 / 10 = %g", + x1), + if(chg0)sprintf(" -> x0 changed too, to x0= %g", x0))) + } + else { ## reasonable R2.rob + if(all(sml.res <- abs(rob.resid) < 4*mad(rob.resid))) { + message("converged in ", it, " iterations") + converged <- TRUE + break + } + ## else + + if (!sml.res[nx]) { ## largest x[] does NOT have small resid + x1 <- exp(max(lx[sml.res])) + if(verbose) cat(sprintf("new x1: log(x1)= %12g\n", log(x1))) + } + else ## throw away small x *ONLY* after large x are fixed + if (!sml.res[1]) { ## smallest x[] does NOT have small resid + x0 <- exp(min(lx[sml.res])) + if(verbose) cat(sprintf("new x0: log(x0)= %12g\n", log(x0))) + } + else ## no change at both ends.. hmm + if(R2.rob > .99999 || it > 10) { + message("** R^2-pseudo-converged in ", it, " iterations") + converged <- NA + break + } + } + if((it <- it+1) > 1000) { + if(verbose) warning("too many iterations!") + break + } + }## end{repeat} + + r <- lm.fit(cbind(1, lx), px) + R2 <- 1 - sum(r$resid^2) / SSY + cf <- r$coefficients + if(do.plot) { + legend("topright", c(sprintf("final int. [%g, %g]", log(x0), log(x1)), + sprintf("coef. = (%g, %g)", cf[1], cf[2])), + pch=NA, bty="n") + abline(v = log(c(x0,x1)), col=adjustcolor(1, 1/2), lty=2) + } + ## browser() + c(cf, df=df, ncp=ncp, converged=converged, + x0=x0, x1=x1, alp=alp, x0.999=x0., R2 = R2, R2.rob=R2.rob) + } > > p.tailAsymp <- function(r, add.central=FALSE, F.add = FALSE) { + x01 <- r[c("x0","x1")] + stopifnot(length(x01) == 2, is.numeric(x01)) + curve(pt(x, df=r["df"], ncp=r["ncp"], log=TRUE,lower.tail=FALSE), + x01[1] / 100, x01[2]*100, log="x", ylab = "", + main = + sprintf(paste("pt(x, df=%g, ncp=%g, log", + if(prod(par("mfrow")> 2)) "..)" + else "=TRUE, lower.tail=FALSE)", sep=''), + r["df"], r["ncp"])) + if(F.add) { + ## This "non-sense" non-central F <==> "two-tail" non-central t + ## i.e., there's NO direct representation *unless* ncp = 0 + colF <- "forestgreen" + curve(pf(x^2, df1=1, df2=r["df"], ncp= r["ncp"]^2, + log=TRUE,lower.tail=FALSE), + add = TRUE, col = colF) + x. <- 10^par("usr")[2] + text(x., pf(x.^2,df1=1, df2=r["df"],ncp= r["ncp"]^2, + log=TRUE,lower.tail=FALSE), + "pf(x^2, (1,df), ncp^2, ...)", adj=c(1,-.2), col = colF) + } + if(add.central) { + curve(pt(x, df=r["df"], log=TRUE,lower.tail=FALSE), add=TRUE, + col = "midnightblue", lwd=2) + } + + abline(v = x01, col ="gray") + + cL <- "tomato" + mtext(sprintf("slope = %-10.6g", r[2]), line=-1 , adj=1, col=cL, cex=par("cex")) + mtext(sprintf("intCpt= %-10.6g", r[1]), line=-2.2, adj=1, col=cL, cex=par("cex")) + abline(a=r[1], b = r[2] * log(10), col = cL) + } > > if(FALSE) + debug(ptRTailAsymp) > r35 <- ptRTailAsymp(df=3, ncp=5) robust R^2 = 0.55639983 'robust' fit not ok ==> x1 := x1 / 10 = 576722 robust R^2 = 0.92602629 'robust' fit not ok ==> x1 := x1 / 10 = 57672.2 robust R^2 = 0.99976074 robust R^2 = 0.99999968 robust R^2 = 0.99999987 robust R^2 = 0.99999988 robust R^2 = 0.99999988 robust R^2 = 0.99999988 robust R^2 = 0.99999995 robust R^2 = 0.99999996 robust R^2 = 0.99999996 robust R^2 = 0.99999996 robust R^2 = 0.99999998 robust R^2 = 0.99999998 robust R^2 = 0.99999998 robust R^2 = 0.99999998 robust R^2 = 0.99999998 robust R^2 = 0.99999999 robust R^2 = 0.99999999 robust R^2 = 0.99999999 robust R^2 = 0.99999999 robust R^2 = 0.99999999 robust R^2 = 0.99999999 robust R^2 = 0.99999999 robust R^2 = 0.99999999 robust R^2 = 0.99999999 robust R^2 = 0.99999999 robust R^2 = 0.99999999 robust R^2 = 0.99999999 robust R^2 = 0.99999999 robust R^2 = 0.99999999 robust R^2 = 0.99999999 robust R^2 = 0.99999999 robust R^2 = 0.99999999 robust R^2 = 0.99999999 robust R^2 = 0.99999999 robust R^2 = 0.99999999 robust R^2 = 0.99999999 robust R^2 = 0.99999999 robust R^2 = 0.99999999 robust R^2 = 0.99999999 robust R^2 = 0.99999999 robust R^2 = 0.99999999 robust R^2 = 0.99999999 robust R^2 = 0.99999999 converged in 45 iterations > if(interactive()) x11(type = "Xlib") > r45 <- ptRTailAsymp(df=4, ncp=5, do.plot=TRUE) robust R^2 = 0.11884457 'robust' fit not ok ==> x1 := x1 / 10 = 349985 robust R^2 = 0.44774985 'robust' fit not ok ==> x1 := x1 / 10 = 34998.5 robust R^2 = 0.92931071 'robust' fit not ok ==> x1 := x1 / 10 = 3499.85 robust R^2 = 0.99999630 new x1: log(x1)= 7.96686 robust R^2 = 0.99999814 new x1: log(x1)= 7.9227 robust R^2 = 0.99999825 new x1: log(x1)= 7.91396 robust R^2 = 0.99999826 new x1: log(x1)= 7.9096 robust R^2 = 0.99999827 new x0: log(x0)= 3.7776 robust R^2 = 0.99999902 new x1: log(x1)= 7.85169 robust R^2 = 0.99999913 new x1: log(x1)= 7.83946 robust R^2 = 0.99999915 new x1: log(x1)= 7.83539 robust R^2 = 0.99999915 new x0: log(x0)= 3.92382 robust R^2 = 0.99999940 new x1: log(x1)= 7.79624 robust R^2 = 0.99999945 new x1: log(x1)= 7.78848 robust R^2 = 0.99999946 new x1: log(x1)= 7.78462 robust R^2 = 0.99999946 new x0: log(x0)= 4.0243 robust R^2 = 0.99999957 new x1: log(x1)= 7.75827 robust R^2 = 0.99999960 new x1: log(x1)= 7.75453 robust R^2 = 0.99999960 new x1: log(x1)= 7.75079 robust R^2 = 0.99999960 new x0: log(x0)= 4.08772 robust R^2 = 0.99999965 new x1: log(x1)= 7.73246 robust R^2 = 0.99999967 new x1: log(x1)= 7.72881 robust R^2 = 0.99999967 new x0: log(x0)= 4.1351 robust R^2 = 0.99999970 new x1: log(x1)= 7.71442 robust R^2 = 0.99999971 new x1: log(x1)= 7.71084 robust R^2 = 0.99999971 new x0: log(x0)= 4.16731 robust R^2 = 0.99999973 new x1: log(x1)= 7.7002 robust R^2 = 0.99999974 new x0: log(x0)= 4.19207 robust R^2 = 0.99999975 new x1: log(x1)= 7.69318 robust R^2 = 0.99999976 new x1: log(x1)= 7.68967 robust R^2 = 0.99999976 new x0: log(x0)= 4.20957 robust R^2 = 0.99999977 new x1: log(x1)= 7.68619 robust R^2 = 0.99999977 new x1: log(x1)= 7.68271 robust R^2 = 0.99999977 new x0: log(x0)= 4.22 robust R^2 = 0.99999978 new x1: log(x1)= 7.67924 robust R^2 = 0.99999978 new x0: log(x0)= 4.23039 robust R^2 = 0.99999978 new x1: log(x1)= 7.67579 robust R^2 = 0.99999979 new x0: log(x0)= 4.23729 robust R^2 = 0.99999979 new x1: log(x1)= 7.67235 robust R^2 = 0.99999979 new x0: log(x0)= 4.24073 robust R^2 = 0.99999979 new x0: log(x0)= 4.24416 robust R^2 = 0.99999979 converged in 42 iterations > r46 <- ptRTailAsymp(df=4, ncp=6, do.plot=TRUE) robust R^2 = 0.30663951 'robust' fit not ok ==> x1 := x1 / 10 = 413611 robust R^2 = 0.46377760 'robust' fit not ok ==> x1 := x1 / 10 = 41361.1 robust R^2 = 0.94373532 'robust' fit not ok ==> x1 := x1 / 10 = 4136.11 robust R^2 = 0.99999792 new x1: log(x1)= 8.19844 robust R^2 = 0.99999845 new x1: log(x1)= 8.16707 robust R^2 = 0.99999849 new x1: log(x1)= 8.16262 robust R^2 = 0.99999850 new x1: log(x1)= 8.15818 robust R^2 = 0.99999851 new x0: log(x0)= 3.971 robust R^2 = 0.99999922 new x1: log(x1)= 8.0995 robust R^2 = 0.99999931 new x1: log(x1)= 8.08297 robust R^2 = 0.99999932 new x0: log(x0)= 4.12741 robust R^2 = 0.99999954 new x1: log(x1)= 8.03941 robust R^2 = 0.99999958 new x1: log(x1)= 8.03158 robust R^2 = 0.99999959 new x1: log(x1)= 8.02767 robust R^2 = 0.99999959 new x1: log(x1)= 8.02377 robust R^2 = 0.99999959 new x0: log(x0)= 4.23661 robust R^2 = 0.99999968 new x1: log(x1)= 7.99723 robust R^2 = 0.99999970 new x1: log(x1)= 7.9897 robust R^2 = 0.99999970 new x0: log(x0)= 4.30799 robust R^2 = 0.99999975 new x1: log(x1)= 7.96759 robust R^2 = 0.99999976 new x1: log(x1)= 7.96393 robust R^2 = 0.99999976 new x0: log(x0)= 4.35923 robust R^2 = 0.99999979 new x1: log(x1)= 7.9495 robust R^2 = 0.99999979 new x1: log(x1)= 7.9459 robust R^2 = 0.99999980 new x0: log(x0)= 4.39154 robust R^2 = 0.99999981 new x1: log(x1)= 7.93523 robust R^2 = 0.99999981 new x1: log(x1)= 7.93168 robust R^2 = 0.99999982 new x0: log(x0)= 4.41989 robust R^2 = 0.99999983 new x1: log(x1)= 7.92465 robust R^2 = 0.99999983 new x1: log(x1)= 7.92114 robust R^2 = 0.99999983 new x0: log(x0)= 4.43741 robust R^2 = 0.99999984 new x1: log(x1)= 7.91765 robust R^2 = 0.99999984 new x1: log(x1)= 7.91417 robust R^2 = 0.99999984 new x0: log(x0)= 4.45134 robust R^2 = 0.99999984 new x0: log(x0)= 4.45827 robust R^2 = 0.99999985 new x1: log(x1)= 7.90725 robust R^2 = 0.99999985 new x0: log(x0)= 4.46517 robust R^2 = 0.99999985 new x0: log(x0)= 4.46862 robust R^2 = 0.99999985 new x1: log(x1)= 7.90381 robust R^2 = 0.99999985 new x0: log(x0)= 4.47206 robust R^2 = 0.99999985 new x0: log(x0)= 4.47549 robust R^2 = 0.99999985 converged in 42 iterations > r410 <- ptRTailAsymp(df=4, ncp=10, do.plot=TRUE) robust R^2 = 0.15002022 'robust' fit not ok ==> x1 := x1 / 10 = 673203 robust R^2 = 0.44951861 'robust' fit not ok ==> x1 := x1 / 10 = 67320.3 robust R^2 = 0.93160095 'robust' fit not ok ==> x1 := x1 / 10 = 6732.03 robust R^2 = 0.99999687 new x1: log(x1)= 8.6118 robust R^2 = 0.99999855 new x1: log(x1)= 8.56333 robust R^2 = 0.99999864 new x1: log(x1)= 8.55025 robust R^2 = 0.99999866 new x1: log(x1)= 8.54591 robust R^2 = 0.99999867 new x0: log(x0)= 4.43518 robust R^2 = 0.99999925 new x1: log(x1)= 8.49241 robust R^2 = 0.99999933 new x1: log(x1)= 8.47617 robust R^2 = 0.99999935 new x0: log(x0)= 4.57676 robust R^2 = 0.99999953 new x1: log(x1)= 8.43714 robust R^2 = 0.99999957 new x1: log(x1)= 8.42941 robust R^2 = 0.99999958 new x1: log(x1)= 8.42555 robust R^2 = 0.99999958 new x0: log(x0)= 4.67307 robust R^2 = 0.99999966 new x1: log(x1)= 8.39926 robust R^2 = 0.99999968 new x1: log(x1)= 8.39553 robust R^2 = 0.99999969 new x1: log(x1)= 8.3918 robust R^2 = 0.99999969 new x0: log(x0)= 4.74008 robust R^2 = 0.99999973 new x1: log(x1)= 8.37352 robust R^2 = 0.99999974 new x1: log(x1)= 8.36989 robust R^2 = 0.99999974 new x1: log(x1)= 8.36625 robust R^2 = 0.99999975 new x0: log(x0)= 4.78364 robust R^2 = 0.99999977 new x1: log(x1)= 8.3555 robust R^2 = 0.99999977 new x1: log(x1)= 8.35192 robust R^2 = 0.99999978 new x0: log(x0)= 4.81935 robust R^2 = 0.99999979 new x1: log(x1)= 8.34131 robust R^2 = 0.99999980 new x0: log(x0)= 4.83698 robust R^2 = 0.99999980 new x1: log(x1)= 8.3378 robust R^2 = 0.99999981 new x1: log(x1)= 8.3343 robust R^2 = 0.99999981 new x0: log(x0)= 4.85449 robust R^2 = 0.99999981 new x1: log(x1)= 8.32733 robust R^2 = 0.99999982 new x0: log(x0)= 4.86492 robust R^2 = 0.99999982 new x1: log(x1)= 8.32387 robust R^2 = 0.99999982 new x0: log(x0)= 4.8753 robust R^2 = 0.99999983 new x1: log(x1)= 8.32041 robust R^2 = 0.99999983 new x0: log(x0)= 4.87875 robust R^2 = 0.99999983 new x0: log(x0)= 4.8822 robust R^2 = 0.99999983 new x1: log(x1)= 8.31697 robust R^2 = 0.99999983 new x0: log(x0)= 4.88563 robust R^2 = 0.99999983 converged in 41 iterations > r420 <- ptRTailAsymp(df=4, ncp=20) robust R^2 = 0.08888495 'robust' fit not ok ==> x1 := x1 / 10 = 1.33222e+06 robust R^2 = 0.44678787 'robust' fit not ok ==> x1 := x1 / 10 = 133222 robust R^2 = 0.92572670 'robust' fit not ok ==> x1 := x1 / 10 = 13322.2 robust R^2 = 0.99999600 robust R^2 = 0.99999860 robust R^2 = 0.99999871 robust R^2 = 0.99999873 robust R^2 = 0.99999874 robust R^2 = 0.99999926 robust R^2 = 0.99999934 robust R^2 = 0.99999935 robust R^2 = 0.99999935 robust R^2 = 0.99999953 robust R^2 = 0.99999957 robust R^2 = 0.99999958 robust R^2 = 0.99999958 robust R^2 = 0.99999966 robust R^2 = 0.99999968 robust R^2 = 0.99999968 robust R^2 = 0.99999968 robust R^2 = 0.99999972 robust R^2 = 0.99999973 robust R^2 = 0.99999973 robust R^2 = 0.99999976 robust R^2 = 0.99999976 robust R^2 = 0.99999977 robust R^2 = 0.99999977 robust R^2 = 0.99999978 robust R^2 = 0.99999979 robust R^2 = 0.99999980 robust R^2 = 0.99999980 robust R^2 = 0.99999980 robust R^2 = 0.99999981 robust R^2 = 0.99999981 robust R^2 = 0.99999982 robust R^2 = 0.99999982 robust R^2 = 0.99999982 robust R^2 = 0.99999982 robust R^2 = 0.99999982 robust R^2 = 0.99999982 robust R^2 = 0.99999983 robust R^2 = 0.99999983 converged in 42 iterations > p.tailAsymp(r420) > p.tailAsymp(r410) > p.tailAsymp(r46) > p.tailAsymp(r45) > > r440 <- ptRTailAsymp(df=4, ncp=40, do.plot=TRUE)## "wrong" robust R^2 = 0.02927432 'robust' fit not ok ==> x1 := x1 / 10 = 4.66178e+06 robust R^2 = 0.15038127 'robust' fit not ok ==> x1 := x1 / 10 = 466178 robust R^2 = 0.36118459 'robust' fit not ok ==> x1 := x1 / 10 = 46617.8 robust R^2 = 0.65582148 'robust' fit not ok ==> x1 := x1 / 10 = 4661.78 robust R^2 = 0.91233348 'robust' fit not ok ==> x1 := x1 / 10 = 466.178 -> x0 changed too, to x0= 233.089 robust R^2 = 0.99342355 converged in 6 iterations > p.tailAsymp(r440) > > r.810 <- ptRTailAsymp(df=.8, ncp=10, do.plot=TRUE) robust R^2 = 0.99980650 new x1: log(x1)= 18.2684 robust R^2 = 0.99999862 new x1: log(x1)= 16.9375 robust R^2 = 0.99999999 new x1: log(x1)= 16.0849 robust R^2 = 1.00000000 new x1: log(x1)= 15.9592 robust R^2 = 1.00000000 new x1: log(x1)= 15.9112 robust R^2 = 1.00000000 ** R^2-pseudo-converged in 6 iterations > p.tailAsymp(r.810)# tip top > > r.518 <- ptRTailAsymp(df=.5, ncp=18, do.plot=TRUE) robust R^2 = 0.99935858 new x0: log(x0)= 18.03 robust R^2 = 0.99989858 new x1: log(x1)= 26.8489 robust R^2 = 0.99982007 new x1: log(x1)= 26.2574 robust R^2 = 0.99982835 new x1: log(x1)= 25.755 robust R^2 = 0.99969431 new x1: log(x1)= 25.4225 robust R^2 = 0.99969607 new x1: log(x1)= 25.1413 robust R^2 = 0.99972195 new x1: log(x1)= 24.8851 robust R^2 = 0.99962719 new x1: log(x1)= 24.7272 robust R^2 = 0.99959265 new x1: log(x1)= 24.573 robust R^2 = 0.99960186 new x1: log(x1)= 24.4289 robust R^2 = 0.99958075 new x1: log(x1)= 24.3201 robust R^2 = 0.99950456 new x1: log(x1)= 24.2571 robust R^2 = 0.99952376 new x1: log(x1)= 24.1823 robust R^2 = 0.99951474 new x1: log(x1)= 24.1145 robust R^2 = 0.99950846 new x1: log(x1)= 24.0475 robust R^2 = 0.99946815 new x1: log(x1)= 24.0114 robust R^2 = 0.99937740 ** R^2-pseudo-converged in 17 iterations > p.tailAsymp(r.518) > ## better (much faster convergence; better range > r.518 <- ptRTailAsymp(df=.5, ncp=18, do.plot=TRUE, f.x0=1/10,f.x1 = 100) robust R^2 = 0.99793674 new x1: log(x1)= 18.3642 robust R^2 = 0.99846461 robust R^2 = 0.99846461 robust R^2 = 0.99846461 robust R^2 = 0.99846461 robust R^2 = 0.99846461 robust R^2 = 0.99846461 robust R^2 = 0.99846461 robust R^2 = 0.99846461 robust R^2 = 0.99846461 robust R^2 = 0.99846461 ** R^2-pseudo-converged in 11 iterations > p.tailAsymp(r.518) > > r.110 <- ptRTailAsymp(df=.1, ncp=10, do.plot=TRUE, f.x0=1e-6,f.x1 = 1e10) robust R^2 = 1.00000000 new x1: log(x1)= 74.517 robust R^2 = 1.00000000 new x1: log(x1)= 74.2879 robust R^2 = 1.00000000 new x1: log(x1)= 74.0838 robust R^2 = 1.00000000 new x1: log(x1)= 73.904 robust R^2 = 1.00000000 new x1: log(x1)= 73.7479 robust R^2 = 1.00000000 new x1: log(x1)= 73.6593 robust R^2 = 1.00000000 new x1: log(x1)= 73.527 robust R^2 = 1.00000000 new x1: log(x1)= 73.3736 robust R^2 = 1.00000000 new x1: log(x1)= 73.2865 robust R^2 = 1.00000000 new x1: log(x1)= 73.1998 robust R^2 = 1.00000000 new x1: log(x1)= 73.0703 robust R^2 = 1.00000000 converged in 12 iterations > p.tailAsymp(r.110) # easy simple > > r71 <- ptRTailAsymp(df=7, ncp=1, do.plot=TRUE) robust R^2 = -0.00346922 'robust' fit not ok ==> x1 := x1 / 10 = 71065.9 robust R^2 = -0.17264703 'robust' fit not ok ==> x1 := x1 / 10 = 7106.59 robust R^2 = 0.45180263 'robust' fit not ok ==> x1 := x1 / 10 = 710.659 robust R^2 = 0.89572260 'robust' fit not ok ==> x1 := x1 / 10 = 71.0659 robust R^2 = 0.99967026 new x0: log(x0)= 2.12006 robust R^2 = 0.99981239 new x0: log(x0)= 2.25095 robust R^2 = 0.99988337 new x0: log(x0)= 2.35974 robust R^2 = 0.99992203 new x0: log(x0)= 2.45122 robust R^2 = 0.99994473 new x0: log(x0)= 2.52923 robust R^2 = 0.99995895 new x0: log(x0)= 2.5952 robust R^2 = 0.99996817 new x0: log(x0)= 2.65198 robust R^2 = 0.99997449 new x0: log(x0)= 2.70038 robust R^2 = 0.99997892 new x0: log(x0)= 2.74263 robust R^2 = 0.99998217 new x0: log(x0)= 2.77917 robust R^2 = 0.99998460 new x0: log(x0)= 2.81186 robust R^2 = 0.99998650 new x0: log(x0)= 2.84092 robust R^2 = 0.99998800 new x0: log(x0)= 2.86656 robust R^2 = 0.99998920 new x0: log(x0)= 2.88893 robust R^2 = 0.99999015 new x0: log(x0)= 2.90957 robust R^2 = 0.99999095 new x0: log(x0)= 2.92719 robust R^2 = 0.99999159 new x0: log(x0)= 2.94325 robust R^2 = 0.99999213 new x0: log(x0)= 2.95779 robust R^2 = 0.99999260 new x0: log(x0)= 2.97086 robust R^2 = 0.99999299 new x0: log(x0)= 2.9825 robust R^2 = 0.99999333 new x0: log(x0)= 2.99404 robust R^2 = 0.99999365 new x0: log(x0)= 3.00421 robust R^2 = 0.99999392 new x0: log(x0)= 3.01304 robust R^2 = 0.99999414 new x0: log(x0)= 3.0218 robust R^2 = 0.99999436 new x0: log(x0)= 3.02926 robust R^2 = 0.99999453 new x0: log(x0)= 3.03667 robust R^2 = 0.99999471 new x0: log(x0)= 3.04281 robust R^2 = 0.99999484 new x0: log(x0)= 3.04892 robust R^2 = 0.99999498 new x0: log(x0)= 3.05378 robust R^2 = 0.99999508 new x0: log(x0)= 3.05863 robust R^2 = 0.99999518 new x0: log(x0)= 3.06345 robust R^2 = 0.99999528 new x0: log(x0)= 3.06826 robust R^2 = 0.99999538 new x0: log(x0)= 3.07185 robust R^2 = 0.99999545 new x0: log(x0)= 3.07543 robust R^2 = 0.99999552 new x0: log(x0)= 3.079 robust R^2 = 0.99999559 new x0: log(x0)= 3.08137 robust R^2 = 0.99999564 new x0: log(x0)= 3.08373 robust R^2 = 0.99999568 new x0: log(x0)= 3.0861 robust R^2 = 0.99999573 new x0: log(x0)= 3.08845 robust R^2 = 0.99999577 new x0: log(x0)= 3.09081 robust R^2 = 0.99999581 new x0: log(x0)= 3.09315 robust R^2 = 0.99999585 new x0: log(x0)= 3.09433 robust R^2 = 0.99999588 new x0: log(x0)= 3.0955 robust R^2 = 0.99999590 new x0: log(x0)= 3.09667 robust R^2 = 0.99999592 new x0: log(x0)= 3.09783 robust R^2 = 0.99999594 new x0: log(x0)= 3.099 robust R^2 = 0.99999596 new x0: log(x0)= 3.10017 robust R^2 = 0.99999598 new x0: log(x0)= 3.10133 robust R^2 = 0.99999600 new x0: log(x0)= 3.10249 robust R^2 = 0.99999602 new x0: log(x0)= 3.10366 robust R^2 = 0.99999604 new x0: log(x0)= 3.10482 robust R^2 = 0.99999606 new x0: log(x0)= 3.10598 robust R^2 = 0.99999608 new x0: log(x0)= 3.10714 robust R^2 = 0.99999610 new x0: log(x0)= 3.10829 robust R^2 = 0.99999612 new x0: log(x0)= 3.10945 robust R^2 = 0.99999614 converged in 60 iterations > p.tailAsymp(r71)# looks good, even though slope = -6.956 > ## --converging to *central* t-distrib: > r7.01 <- ptRTailAsymp(df=7, ncp=.01, do.plot=TRUE) robust R^2 = -0.05986955 'robust' fit not ok ==> x1 := x1 / 10 = 48067.2 robust R^2 = -0.36332463 'robust' fit not ok ==> x1 := x1 / 10 = 4806.72 robust R^2 = 0.54057131 'robust' fit not ok ==> x1 := x1 / 10 = 480.672 robust R^2 = 0.97553626 'robust' fit not ok ==> x1 := x1 / 10 = 48.0672 robust R^2 = 0.99928907 new x0: log(x0)= 1.71522 robust R^2 = 0.99956467 new x0: log(x0)= 1.83832 robust R^2 = 0.99971650 new x0: log(x0)= 1.94217 robust R^2 = 0.99980432 new x0: log(x0)= 2.02913 robust R^2 = 0.99985740 new x0: log(x0)= 2.10478 robust R^2 = 0.99989219 new x0: log(x0)= 2.16849 robust R^2 = 0.99991508 new x0: log(x0)= 2.22308 robust R^2 = 0.99993095 new x0: log(x0)= 2.27096 robust R^2 = 0.99994251 new x0: log(x0)= 2.31264 robust R^2 = 0.99995106 new x0: log(x0)= 2.34856 robust R^2 = 0.99995744 new x0: log(x0)= 2.3806 robust R^2 = 0.99996246 new x0: log(x0)= 2.40897 robust R^2 = 0.99996643 new x0: log(x0)= 2.43388 robust R^2 = 0.99996959 new x0: log(x0)= 2.45548 robust R^2 = 0.99997210 new x0: log(x0)= 2.47534 robust R^2 = 0.99997423 new x0: log(x0)= 2.49212 robust R^2 = 0.99997591 new x0: log(x0)= 2.50733 robust R^2 = 0.99997734 new x0: log(x0)= 2.52099 robust R^2 = 0.99997856 new x0: log(x0)= 2.53317 robust R^2 = 0.99997959 new x0: log(x0)= 2.54389 robust R^2 = 0.99998046 new x0: log(x0)= 2.5532 robust R^2 = 0.99998119 new x0: log(x0)= 2.56113 robust R^2 = 0.99998179 new x0: log(x0)= 2.56901 robust R^2 = 0.99998237 new x0: log(x0)= 2.57553 robust R^2 = 0.99998283 new x0: log(x0)= 2.58202 robust R^2 = 0.99998329 new x0: log(x0)= 2.58719 robust R^2 = 0.99998364 new x0: log(x0)= 2.59234 robust R^2 = 0.99998398 new x0: log(x0)= 2.59746 robust R^2 = 0.99998432 new x0: log(x0)= 2.60129 robust R^2 = 0.99998456 new x0: log(x0)= 2.60511 robust R^2 = 0.99998480 new x0: log(x0)= 2.60892 robust R^2 = 0.99998504 new x0: log(x0)= 2.61145 robust R^2 = 0.99998520 new x0: log(x0)= 2.61397 robust R^2 = 0.99998535 new x0: log(x0)= 2.61649 robust R^2 = 0.99998550 new x0: log(x0)= 2.61901 robust R^2 = 0.99998565 new x0: log(x0)= 2.62026 robust R^2 = 0.99998573 new x0: log(x0)= 2.62151 robust R^2 = 0.99998580 new x0: log(x0)= 2.62277 robust R^2 = 0.99998588 new x0: log(x0)= 2.62402 robust R^2 = 0.99998595 new x0: log(x0)= 2.62527 robust R^2 = 0.99998602 new x0: log(x0)= 2.62652 robust R^2 = 0.99998609 new x0: log(x0)= 2.62776 robust R^2 = 0.99998616 new x0: log(x0)= 2.62901 robust R^2 = 0.99998624 new x0: log(x0)= 2.63025 robust R^2 = 0.99998631 new x0: log(x0)= 2.6315 robust R^2 = 0.99998638 new x0: log(x0)= 2.63274 robust R^2 = 0.99998645 new x0: log(x0)= 2.63398 robust R^2 = 0.99998652 new x0: log(x0)= 2.63522 robust R^2 = 0.99998659 converged in 53 iterations > p.tailAsymp(r7.01)# looks ok > r7.001 <- ptRTailAsymp(df=7, ncp=.001, do.plot=TRUE) robust R^2 = -0.00776460 'robust' fit not ok ==> x1 := x1 / 10 = 47874.3 robust R^2 = -0.26070239 'robust' fit not ok ==> x1 := x1 / 10 = 4787.43 robust R^2 = 0.80057238 'robust' fit not ok ==> x1 := x1 / 10 = 478.743 robust R^2 = 0.94533935 'robust' fit not ok ==> x1 := x1 / 10 = 47.8743 robust R^2 = 0.99928313 new x0: log(x0)= 1.7112 robust R^2 = 0.99956091 new x0: log(x0)= 1.8343 robust R^2 = 0.99971401 new x0: log(x0)= 1.93815 robust R^2 = 0.99980258 new x0: log(x0)= 2.0251 robust R^2 = 0.99985612 new x0: log(x0)= 2.10076 robust R^2 = 0.99989122 new x0: log(x0)= 2.16447 robust R^2 = 0.99991432 new x0: log(x0)= 2.21905 robust R^2 = 0.99993034 new x0: log(x0)= 2.26694 robust R^2 = 0.99994201 new x0: log(x0)= 2.30862 robust R^2 = 0.99995063 new x0: log(x0)= 2.34454 robust R^2 = 0.99995707 new x0: log(x0)= 2.37657 robust R^2 = 0.99996214 new x0: log(x0)= 2.40495 robust R^2 = 0.99996615 new x0: log(x0)= 2.42986 robust R^2 = 0.99996934 new x0: log(x0)= 2.45146 robust R^2 = 0.99997187 new x0: log(x0)= 2.47132 robust R^2 = 0.99997402 new x0: log(x0)= 2.4895 robust R^2 = 0.99997585 new x0: log(x0)= 2.50469 robust R^2 = 0.99997729 new x0: log(x0)= 2.51834 robust R^2 = 0.99997851 new x0: log(x0)= 2.5305 robust R^2 = 0.99997955 new x0: log(x0)= 2.54122 robust R^2 = 0.99998042 new x0: log(x0)= 2.55052 robust R^2 = 0.99998115 new x0: log(x0)= 2.55976 robust R^2 = 0.99998185 new x0: log(x0)= 2.56762 robust R^2 = 0.99998243 new x0: log(x0)= 2.57543 robust R^2 = 0.99998298 new x0: log(x0)= 2.5819 robust R^2 = 0.99998343 new x0: log(x0)= 2.58834 robust R^2 = 0.99998387 new x0: log(x0)= 2.59347 robust R^2 = 0.99998420 new x0: log(x0)= 2.59857 robust R^2 = 0.99998453 new x0: log(x0)= 2.60239 robust R^2 = 0.99998478 new x0: log(x0)= 2.60619 robust R^2 = 0.99998502 new x0: log(x0)= 2.60998 robust R^2 = 0.99998525 new x0: log(x0)= 2.6125 robust R^2 = 0.99998540 new x0: log(x0)= 2.61501 robust R^2 = 0.99998556 new x0: log(x0)= 2.61752 robust R^2 = 0.99998570 new x0: log(x0)= 2.62003 robust R^2 = 0.99998585 new x0: log(x0)= 2.62253 robust R^2 = 0.99998600 new x0: log(x0)= 2.62502 robust R^2 = 0.99998615 new x0: log(x0)= 2.62627 robust R^2 = 0.99998622 new x0: log(x0)= 2.62751 robust R^2 = 0.99998629 new x0: log(x0)= 2.62875 robust R^2 = 0.99998636 new x0: log(x0)= 2.63 robust R^2 = 0.99998643 new x0: log(x0)= 2.63123 robust R^2 = 0.99998650 new x0: log(x0)= 2.63247 robust R^2 = 0.99998657 new x0: log(x0)= 2.63371 robust R^2 = 0.99998664 new x0: log(x0)= 2.63495 robust R^2 = 0.99998671 new x0: log(x0)= 2.63618 robust R^2 = 0.99998678 new x0: log(x0)= 2.63742 robust R^2 = 0.99998684 new x0: log(x0)= 2.63865 robust R^2 = 0.99998691 converged in 53 iterations > p.tailAsymp(r7.001, add.central=TRUE)# --- but the slope is not quite -7 ... > ## hmm: > tt <- lseq(10,100, length=1000) > px <- pt(tt, df=7, log=TRUE,lower.tail=FALSE) > plot(px ~ log(tt), pch=".", main = "pt(tt, df=7, log=TRUE, lower.t=F)") > ll <- lm(px ~ log(tt)) > summary(ll) Call: lm(formula = px ~ log(tt)) Residuals: Min 1Q Median 3Q Max -0.084028 -0.017043 0.006405 0.022235 0.027875 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 4.595546 0.004271 1076 <2e-16 *** log(tt) -6.930061 0.001214 -5707 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.02555 on 998 degrees of freedom Multiple R-squared: 1, Adjusted R-squared: 1 F-statistic: 3.257e+07 on 1 and 998 DF, p-value: < 2.2e-16 > ## Estimate Std. Error t value Pr(>|t|) > ## (Intercept) 4.595546 0.004271 1076 <2e-16 *** > ## log(tt) -6.930061 0.001214 -5707 <2e-16 *** > ## ^^^^^^^^^ ....... > ##--- significantly different from -7 > > tt <- lseq(100,1e7, length=1000) # much larger range: t -> oo > px <- pt(tt, df=7, log=TRUE,lower.tail=FALSE) > plot(px ~ log(tt), pch=".", + main = "pt(tt, df=7, log=TRUE, lower.t=F) -- larger range(tt)") > ll <- lm(px ~ log(tt)) > summary(ll)##--> clearly: slope == -7 -- ok Call: lm(formula = px ~ log(tt)) Residuals: Min 1Q Median 3Q Max -1.820e-03 -7.448e-05 4.739e-05 1.620e-04 2.307e-04 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 4.883e+00 2.802e-05 174231 <2e-16 *** log(tt) -7.000e+00 2.575e-06 -2718316 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.0002709 on 998 degrees of freedom Multiple R-squared: 1, Adjusted R-squared: 1 F-statistic: 7.389e+12 on 1 and 998 DF, p-value: < 2.2e-16 > > ## Left tail > tt <- rev(- lseq(100,1e7, length=1000)) > px <- pt(tt, df=4.5, log=TRUE) > plot(px ~ log(-tt), pch=".", + main = "pt(tt, df=7, log=TRUE, lower.tail=TRUE)") > ll <- lm(px ~ log(-tt)) > summary(ll)##--> clearly: slope == -4.5 -- ok Call: lm(formula = px ~ log(-tt)) Residuals: Min 1Q Median 3Q Max -7.160e-04 -2.930e-05 1.864e-05 6.375e-05 9.078e-05 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.658e+00 1.103e-05 150365 <2e-16 *** log(-tt) -4.500e+00 1.013e-06 -4441649 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.0001066 on 998 degrees of freedom Multiple R-squared: 1, Adjusted R-squared: 1 F-statistic: 1.973e+13 on 1 and 998 DF, p-value: < 2.2e-16 > > ### ====> (empirical) Theorem: slope == -df > ### ============ > > ### and same is true for left tail ---- see below > (rmat0 <- t(sapply(ls(patt="^r.*[0-9]$"), get))) lx df ncp converged x0 x1 r.110 -0.3632572 -0.1000000 0.1 10.000 1 2.644807e+22 5.420146e+31 r.518 1.2460979 -0.5032407 0.5 18.000 NA 1.095266e+06 9.450517e+07 r.810 1.5944075 -0.8000039 0.8 10.000 NA 4.125975e+04 8.130854e+06 r35 5.2625625 -2.9996318 3.0 5.000 1 1.452059e+02 7.465496e+03 r410 9.9450196 -3.9975038 4.0 10.000 1 1.323744e+02 4.092752e+03 r420 12.6719920 -3.9974572 4.0 20.000 1 2.546093e+02 7.814041e+03 r440 0.7737503 -0.9166890 4.0 40.000 1 2.330892e+02 4.661784e+02 r45 7.3326766 -3.9972146 4.0 5.000 1 6.969736e+01 2.148118e+03 r46 8.0013752 -3.9976673 4.0 6.000 1 8.783782e+01 2.707579e+03 r7.001 4.6028639 -6.9255854 7.0 0.001 1 1.399427e+01 4.787431e+01 r7.01 4.6265655 -6.9253751 7.0 0.010 1 1.394639e+01 4.806725e+01 r71 7.1954054 -6.9558651 7.0 1.000 1 2.240873e+01 7.106591e+01 alp x0.999 R2 R2.rob r.110 0.999 2.644807e+28 1.0000000 1.0000000 r.518 0.999 1.095266e+07 0.9992306 0.9984646 r.810 0.999 4.125975e+04 1.0000000 1.0000000 r35 0.999 5.767216e+01 1.0000000 1.0000000 r410 0.999 6.732026e+01 0.9999999 0.9999998 r420 0.999 1.332224e+02 0.9999999 0.9999998 r440 0.992 4.661784e+02 0.9938396 0.9934236 r45 0.999 3.499847e+01 0.9999999 0.9999998 r46 0.999 4.136106e+01 0.9999999 0.9999999 r7.001 0.999 4.787431e+00 0.9999892 0.9999869 r7.01 0.999 4.806725e+00 0.9999889 0.9999866 r71 0.999 7.106591e+00 0.9999969 0.9999961 > rmat <- rmat0[ rmat0[,"R2.rob"] > 0.9999 ,] > rmat[,2:3] lx df r.110 -0.1000000 0.1 r.810 -0.8000039 0.8 r35 -2.9996318 3.0 r410 -3.9975038 4.0 r420 -3.9974572 4.0 r45 -3.9972146 4.0 r46 -3.9976673 4.0 r7.001 -6.9255854 7.0 r7.01 -6.9253751 7.0 r71 -6.9558651 7.0 > > p.tailAsymp(r2.30 <- ptRTailAsymp(df= 2, ncp=30)) # very nice robust R^2 = 0.99497006 robust R^2 = 0.99999989 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 robust R^2 = 1.00000000 converged in 41 iterations > ## "Large" ncp : (the search and resulting slope seem *wrong* ..) > p.tailAsymp(r2.50 <- ptRTailAsymp(df= 2, ncp=50, do.plot=TRUE)) robust R^2 = 0.03003512 'robust' fit not ok ==> x1 := x1 / 10 = 4.38697e+06 robust R^2 = 0.15184065 'robust' fit not ok ==> x1 := x1 / 10 = 438697 robust R^2 = 0.36358778 'robust' fit not ok ==> x1 := x1 / 10 = 43869.7 robust R^2 = 0.65867858 'robust' fit not ok ==> x1 := x1 / 10 = 4386.97 robust R^2 = 0.91379634 'robust' fit not ok ==> x1 := x1 / 10 = 438.697 -> x0 changed too, to x0= 219.348 robust R^2 = 0.99384916 converged in 6 iterations > > showProc.time() Time (user system elapsed): 5.64 0.28 5.94 > > ## Now do a larger set: > > nd <- length(df. <- c(.1, .8, 1, 1.5, 2, 4, 10)) > nc <- length(nc. <- c(.01, .1, 1, 2, 5)) > ## for indexing etc: > c.df <- formatC(df., width=1) > c.nc <- formatC(nc., width=1) > c.c <- outer(c.df, c.nc, paste, sep="_") > indR <- function(i.df, i.nc) paste(c.df[i.df], c.nc[i.nc], sep="_") > > sfil1 <- file.path(sdir, "pnt-prec-sim1.rds") > if(!doExtras && file.exists(sfil1)) { + + Res <- readRDS_(sfil1) + + } else { ## do run the simulation [always if(doExtras)] : --------------- + + r35 <- ptRTailAsymp(df=3, ncp=5) + + if(names(r35)[1] == "") names(r35)[1] <- "intercpt" + Res <- matrix(NA_real_, nrow = length(r35), ## <-- length of output + ncol = nd*nc, + dimnames=list(names(r35), c.c)) + for(i.df in seq_along(df.)) { + df <- df.[i.df] + cat("\ndf = ", formatC(df)," : \n----------\n") + for(i.nc in seq_along(nc.)) { + cat(i.nc,"") + ncp <- nc.[i.nc] + r <- try(ptRTailAsymp(df=df, ncp=ncp)) + Res[, indR(i.df,i.nc) ] <- if(inherits(r, "try-error")) NA else r + } + } + attr(Res, "version") <- list(R.version) + save2RDS(Res, file=sfil1) + }##-- end{ run simulation } -------------------------------------------- Reading from D:/RCompile/CRANincoming/R-devel/lib/DPQ/safe/pnt-prec-sim1.rds Time (user system elapsed): 0 0 0 > > dR <- as.data.frame(t(Res)) > > op <- mult.fig(mfrow=c(nd,nc),marP=-1)$old.par > op2 <- par(cex = par("cex") * 3/4) # or even more > for(i.df in seq_along(df.)) + for(i.nc in seq_along(nc.)) { + r <- Res[, indR(i.df,i.nc) ] + if(is.na(r["df"])) { frame() } else p.tailAsymp(r) + # ------------ + } > par(op); par(op2) > showProc.time() Time (user system elapsed): 0.08 0 0.07 > > > ### ====> Theorem: slope == -df > ### ============ [now have prove for ncp=0: central t] > ### and same is true for left tail > ## but I see extra discontinuous behavior at x = 0 : > > ## df = 10 and various ncp : > curve(pt(x, df=10, ncp=10, log=TRUE), -250,50, ylim=c(-50,0), n=10001) > title("pt(x, df=10, ncp= *, log=TRUE), -250,50, ..) for various ncp") > ## > curve(pt(x, df=10, ncp=0.01, log=TRUE), add=TRUE, n=10001, col=2) There were 816 warnings (use warnings() to see them) > curve(pt(x, df=10, ncp=0.0, log=TRUE), add=TRUE, n=10001, col=2) > ## > curve(pt(x, df=10, ncp=0.1, log=TRUE), add=TRUE, n=10001, col="red2") There were 790 warnings (use warnings() to see them) > curve(pt(x, df=10, ncp=0.2, log=TRUE), add=TRUE, n=10001, col="red2") There were 762 warnings (use warnings() to see them) > curve(pt(x, df=10, ncp=0.5, log=TRUE), add=TRUE, n=10001, col="red2") There were 674 warnings (use warnings() to see them) > ## > curve(pt(x, df=10, ncp=1, log=TRUE), add=TRUE, n=10001, col="blue2") There were 519 warnings (use warnings() to see them) > curve(pt(x, df=10, ncp=2, log=TRUE), add=TRUE, n=10001, col="blue2") There were 182 warnings (use warnings() to see them) > curve(pt(x, df=10, ncp=5, log=TRUE), add=TRUE, n=10001, col="blue2") > > > ## same, df=5: here we see increasing problem at x=0 for ncp -> large > curve(pt(x, df=5, ncp=100, log=TRUE), -5000,300, ylim=c(-50,0), n=10001) > curve(pt(x, df=5, ncp=0, log=TRUE), add=TRUE, n=10001, + col="blue",lty=3, lwd=3) > for(n. in c(1:8) / 10) + curve(pt(x, df=5, ncp=n., log=TRUE), add=TRUE, n=2001, col="blue2") There were 330 warnings (use warnings() to see them) > for(n. in c(1:8)) + curve(pt(x, df=5, ncp=n., log=TRUE), add=TRUE, n=2001, col="red2") There were 23 warnings (use warnings() to see them) > for(n. in 10*c(1:8)) + curve(pt(x, df=5, ncp=n., log=TRUE), add=TRUE, n=2001, col="green2") > > showProc.time() Time (user system elapsed): 0.83 0 0.83 > > 1 [1] 1 > ## From: Michael E Meredith > ## To: Martin Maechler > ## Subject: Re: [Rd] Noncentral dt() with tiny 'x' values (PR#8874) > ## Date: Thu, 18 May 2006 19:54:09 +0800 > > ## Hello Martin, > > ## Thanks for the response. I see what you mean about the plots! > > ## One thing I should also mention is that I'm getting huge numbers of warnings: > > ## " full precision was not achieved in 'pnt' " > > ## I've been hunting for a pbm in my code, but have now found I get it > ## with 'power.t.test(*, sd = NULL) > # eg. > (ptt <- power.t.test(20, 1, power=0.8, sd=NULL, tol = 1e-8)) Two-sample t test power calculation n = 20 delta = 1 sd = 1.099953 sig.level = 0.05 power = 0.8 alternative = two.sided NOTE: n is number in *each* group > stopifnot(inherits(ptt, "power.htest"), is.list(ptt), + all.equal(1.0999525, ptt$sd, tol = 1e-7)) > ## MM: This has been fixed for R 2.4.0, with NEWS entry > ## > ## o pt() with a very small (or zero) non-centrality parameter could > ## give an unduly stringent warning about 'full precision was not > ## achieved'. (PR#9171) > > > > ## I presume this is another related old buglet which is now appearing > ## due to the improved reporting of warnings in R 2.3.0. > > ## Regards, Mike. > > > ## MM: The note below is > ## > ## Fixed for R 2.3.1 --- with NEWS entry > ## ----- -------- > ## o dt(x, df, ncp= not.0) no longer gives erratic values for > ## |x| < ~1e-12. (PR#8874) > > ## At 17:57 18/05/2006, you wrote: > ## > >>>>> "MikeMer" == meredith > ## > >>>>> on Thu, 18 May 2006 03:52:51 +0200 (CEST) writes: > ## > > ## > MikeMer> Full_Name: Mike Meredith > ## > MikeMer> Version: 2.3.0 > ## > MikeMer> OS: WinXP SP2 > ## > MikeMer> Submission from: (NULL) (210.195.228.29) > ## > > ## > > ## > > ## > MikeMer> Using dt() with a non-centrality parameter and > ## > near-zero values for 'x' results > ## > MikeMer> in erratic output. Try this: > > tst <- c(1e-12, 1e-13, 1e-14, 1e-15, 1e-16, 1e-17, 0) > dt(tst,16,1) [1] 0.2382217 0.2382217 0.2382217 0.2382217 0.2382217 0.2382217 0.2382217 > > ## > MikeMer> I get: 0.2381019 0.2385462 0.2296557 0.1851817 > ## > 0.6288373 3.8163916 (!!) > ## > MikeMer> 0.2382217 > ## > > ## >I get quite different values (several '0', BTW), which just > ## >confirms the erratic nature. > ## > > ## >As often, plots give even a clearer picture: > > > x <- lseq(1e-3, 1e-33, length= 301) > > plot(x, dt(x, df=16, ncp=1), type = "o", cex=.5, log = "x") > plot(x, dt(x, df=16, ncp=0.1), type = "o", cex=.5, log = "x") > plot(x, dt(x, df= 3, ncp=0.1), type = "o", cex=.5, log = "x") > > ## > > ## > > ## > MikeMer> The 0.238 values are okay, the others nonsense, and > ## > MikeMer> they cause confusing spikes on plots of dt() vs 'x' > ## > MikeMer> if 'x' happens to include tiny values. (Other > ## > MikeMer> values of df and ncp also malfunction, but not all > ## > MikeMer> give results out by an order of magnitude!) > ## > > ## >I think almost all do, once you start looking at plots like the above. > ## > > ## > MikeMer> I'm using the work-around dt(round(x,10),...), but > ## > MikeMer> dt() should really take care of this itself. > ## > > ## >or actually rather do something more smart; the cutoff at 1e-10 > ## >is quite crude. > ## > > ## >Note that this is not a new bug at all; but rather as old as > ## >we have dt(*, ncp= .) in R. > ## > > ## >Thanks for reporting it! > ## >Martin Maechler, ETH Zurich > > ## =============================================== > ## Michael E Meredith > ## Wildlife Conservation Society (WCS) Malaysia Program > ## 7 Jalan Ridgeway, 93250 Kuching, Malaysia > ## Fax: +60-82-252 799 Mobile: +60-19 888 1533 > ## email: meredith@easynet.co.uk http://www.mered.org.uk > ## Program website: http://www.wcsmalaysia.org > > > proc.time() user system elapsed 7.17 0.45 7.62