
R version 3.6.2 Patched (2020-02-12 r77795) -- "Dark and Stormy Night"
Copyright (C) 2020 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

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.

> ####	d|ensity
> ####	p|robability (cumulative)
> ####	q|uantile
> ####	r|andom number generation
> ####
> ####	Functions for  ``d/p/q/r''
> 
> .ptime <- proc.time()
> F <- FALSE
> T <- TRUE
> 
> options(warn = 2)
> ##      ======== No warnings, unless explicitly asserted via
> assertWarning <- tools::assertWarning
> 
> as.nan <- function(x) { x[is.na(x) & !is.nan(x)] <- NaN ; x }
> ###-- these are identical in ./arith-true.R ["fixme": use source(..)]
> opt.conformance <- 0
> Meps <- .Machine $ double.eps
> xMax <- .Machine $ double.xmax
> options(rErr.eps = 1e-30)
> rErr <- function(approx, true, eps = getOption("rErr.eps", 1e-30))
+ {
+     ifelse(Mod(true) >= eps,
+ 	   1 - approx / true, # relative error
+ 	   true - approx)     # absolute error (e.g. when true=0)
+ }
> ## Numerical equality: Here want "rel.error" almost always:
> All.eq <- function(x,y) {
+     all.equal.numeric(x,y, tolerance = 64*.Machine$double.eps,
+                       scale = max(0, mean(abs(x), na.rm=TRUE)))
+ }
> if(!interactive())
+     set.seed(123)
> 
> ## The prefixes of ALL the PDQ & R functions
> PDQRinteg <- c("binom", "geom", "hyper", "nbinom", "pois","signrank","wilcox")
> PDQR <- c(PDQRinteg, "beta", "cauchy", "chisq", "exp", "f", "gamma",
+ 	  "lnorm", "logis", "norm", "t","unif","weibull")
> PQonly <- c("tukey")
> 
> ###--- Discrete Distributions --- Consistency Checks  pZZ = cumsum(dZZ)
> 
> ##for(pre in PDQRinteg) { n <- paste("d",pre,sep=""); cat(n,": "); str(get(n))}
> 
> ##__ 1. Binomial __
> 
> ## Cumulative Binomial '==' Cumulative F :
> ## Abramowitz & Stegun, p.945-6;  26.5.24  AND	26.5.28 :
> n0 <- 50; n1 <- 16; n2 <- 20; n3 <- 8
> for(n in rbinom(n1, size = 2*n0, p = .4)) {
+     for(p in c(0,1,rbeta(n2, 2,4))) {
+ 	for(k in rbinom(n3, size = n,  prob = runif(1)))
+ 	    ## For X ~ Bin(n,p), compute 1 - P[X > k] = P[X <= k] in three ways:
+ 	    stopifnot(all.equal(       pbinom(0:k, size = n, prob = p),
+ 				cumsum(dbinom(0:k, size = n, prob = p))),
+ 		      all.equal(if(k==n || p==0) 1 else
+ 				pf((k+1)/(n-k)*(1-p)/p, df1=2*(n-k), df2=2*(k+1)),
+ 				sum(dbinom(0:k, size = n, prob = p))))
+     }
+ }
> 
> ##__ 2. Geometric __
> for(pr in seq(1e-10,1,len=15)) # p=0 is not a distribution
+     stopifnot(All.eq((dg <- dgeom(0:10, pr)),
+ 		     pr * (1-pr)^(0:10)),
+ 	      All.eq(cumsum(dg), pgeom(0:10, pr)))
> 
> 
> ##__ 3. Hypergeometric __
> 
> m <- 10; n <- 7
> for(k in 2:m) {
+     x <- 0:(k+1)
+     stopifnot(All.eq(phyper(x, m, n, k), cumsum(dhyper(x, m, n, k))))
+ }
> 
> ##__ 4. Negative Binomial __
> 
> ## PR #842
> for(size in seq(0.8,2, by=.1))
+     stopifnot(all.equal(cumsum(dnbinom(0:7, size, .5)),
+ 			       pnbinom(0:7, size, .5)))
> stopifnot(All.eq(pnbinom(c(1,3), .9, .5),
+ 		 c(0.777035760338812, 0.946945347071519)))
> 
> ##__ 5. Poisson __
> 
> stopifnot(dpois(0:5,0)		 == c(1, rep(0,5)),
+ 	  dpois(0:5,0, log=TRUE) == c(0, rep(-Inf, 5)))
> 
> ## Cumulative Poisson '==' Cumulative Chi^2 :
> ## Abramowitz & Stegun, p.941 :	 26.4.21 (26.4.2)
> n1 <- 20; n2 <- 16
> for(lambda in rexp(n1))
+     for(k in rpois(n2, lambda))
+ 	stopifnot(all.equal(1 - pchisq(2*lambda, 2*(1+ 0:k)),
+ 			    pp <- cumsum(dpois(0:k, lambda=lambda)),
+ 			    tolerance = 100*Meps),
+ 		  all.equal(pp, ppois(0:k, lambda=lambda), tolerance = 100*Meps),
+ 		  all.equal(1 - pp, ppois(0:k, lambda=lambda, lower.tail = FALSE)))
> 
> 
> ##__ 6. SignRank __
> for(n in rpois(32, lam=8)) {
+     x <- -1:(n + 4)
+     stopifnot(All.eq(psignrank(x, n), cumsum(dsignrank(x, n))))
+ }
> 
> ##__ 7. Wilcoxon (symmetry & cumulative) __
> is.sym <- TRUE
> for(n in rpois(5, lam=6))
+     for(m in rpois(15, lam=8)) {
+ 	x <- -1:(n*m + 1)
+ 	fx <- dwilcox(x, n, m)
+ 	Fx <- pwilcox(x, n, m)
+ 	is.sym <- is.sym & all(fx == dwilcox(x, m, n))
+ 	stopifnot(All.eq(Fx, cumsum(fx)))
+     }
> stopifnot(is.sym)
> 
> 
> ###-------- Continuous Distributions ----------
> 
> ##---  Gamma (incl. central chi^2) Density :
> x <- round(rgamma(100, shape = 2),2)
> for(sh in round(rlnorm(30),2)) {
+     Ga <- gamma(sh)
+     for(sig in round(rlnorm(30),2))
+ 	stopifnot(all.equal((d1 <- dgamma(	 x,   shape = sh, scale = sig)),
+                             (d2 <- dgamma(x/sig, shape = sh, scale = 1) / sig),
+                             tolerance = 1e-14)## __ad interim__ was 1e-15
+                   ,
+                   All.eq(d1, (d3 <- 1/(Ga * sig^sh) * x^(sh-1) * exp(-x/sig)))
+                   )
+ }
> 
> stopifnot(pgamma(1,Inf,scale=Inf) == 0)
> ## Also pgamma(Inf,Inf) == 1 for which NaN was slightly more appropriate
> assertWarning(stopifnot(
+     is.nan(c(pgamma(Inf,  1,scale=Inf),
+              pgamma(Inf,Inf,scale=Inf)))))
> scLrg <- c(2,100, 1e300*c(.1, 1,10,100), 1e307, xMax, Inf)
> stopifnot(pgamma(Inf, 1, scale=xMax) == 1,
+           pgamma(xMax,1, scale=Inf) == 0,
+           all.equal(pgamma(1e300, 2, scale= scLrg, log=TRUE),
+                     c(0, 0, -0.000499523968713701, -1.33089326820406,
+                       -5.36470502873211, -9.91015144019122,
+                       -32.9293385491433, -38.707517174609, -Inf),
+                     tolerance = 2e-15)
+           )
> 
> p <- 7e-4; df <- 0.9
> stopifnot(
+ abs(1-c(pchisq(qchisq(p, df),df)/p, # was 2.31e-8 for R <= 1.8.1
+         pchisq(qchisq(1-p, df,lower=FALSE),df,lower=FALSE)/(1-p),# was 1.618e-11
+         pchisq(qchisq(log(p), df,log=TRUE),df, log=TRUE)/log(p), # was 3.181e-9
+         pchisq(qchisq(log1p(-p),df,log=T,lower=F),df, log=T,lower=F)/log1p(-p)
+         )# 32b-i386: (2.2e-16, 0,0, 3.3e-16); Opteron: (2.2e-16, 0,0, 2.2e-15)
+     ) < 1e-14
+ )
> 
> ##-- non central Chi^2 :
> xB <- c(2000,1e6,1e50,Inf)
> for(df in c(0.1, 1, 10))
+     for(ncp in c(0, 1, 10, 100)) stopifnot(pchisq(xB, df=df, ncp=ncp) == 1)
> stopifnot(all.equal(qchisq(0.025,31,ncp=1,lower.tail=FALSE),# inf.loop PR#875
+                     49.7766246561514, tolerance = 1e-11))
> for(df in c(0.1, 0.5, 1.5, 4.7, 10, 20,50,100)) {
+     xx <- c(10^-(5:1), .9, 1.2, df + c(3,7,20,30,35,38))
+     pp <- pchisq(xx, df=df, ncp = 1) #print(pp)
+     dtol <- 1e-12 *(if(2 < df && df <= 50) 64 else if(df > 50) 20000 else 501)
+     stopifnot(all.equal(xx, qchisq(pp, df=df, ncp=1), tolerance = dtol))
+ }
> 
> ## p ~= 1 (<==> 1-p ~= 0) -- gave infinite loop in R <= 1.8.1 -- PR#6421
> psml <- 2^-(10:54)
> q0 <- qchisq(psml,    df=1.2, ncp=10, lower.tail=FALSE)
> q1 <- qchisq(1-psml, df=1.2, ncp=10) # inaccurate in the tail
> p0 <- pchisq(q0, df=1.2, ncp=10, lower.tail=FALSE)
> p1 <- pchisq(q1, df=1.2, ncp=10, lower.tail=FALSE)
> iO <- 1:30
> stopifnot(all.equal(q0[iO], q1[iO], tolerance = 1e-5),# 9.86e-8
+           all.equal(p0[iO], psml[iO])) # 1.07e-13
> 
> ##--- Beta (need more):
> 
> ## big a & b (PR #643)
> stopifnot(is.finite(a <- rlnorm(20, 5.5)), a > 0,
+           is.finite(b <- rlnorm(20, 6.5)), b > 0)
> pab <- expand.grid(seq(0,1,by=.1), a, b)
> p <- pab[,1]; a <- pab[,2]; b <- pab[,3]
> stopifnot(all.equal(dbeta(p,a,b),
+                     exp(pab <- dbeta(p,a,b, log = TRUE)), tolerance = 1e-11))
> sp <- sample(pab, 50)
> if(!interactive())
+ stopifnot(which(isI <- sp == -Inf) ==
+               c(3, 10, 14, 18, 24, 32, 35, 41, 42, 45, 46, 47),
+           all.equal(range(sp[!isI]), c(-2888.393250, 3.181137))
+           )
> 
> 
> ##--- Normal (& Lognormal) :
> 
> stopifnot(
+     qnorm(0) == -Inf, qnorm(-Inf, log = TRUE) == -Inf,
+     qnorm(1) ==  Inf, qnorm( 0,   log = TRUE) ==  Inf)
> 
> assertWarning(stopifnot(
+     is.nan(qnorm(1.1)),
+     is.nan(qnorm(-.1))))
> 
> x <- c(-Inf, -1e100, 1:6, 1e200, Inf)
> stopifnot(
+     dnorm(x,3,s=0) == c(0,0,0,0, Inf, 0,0,0,0,0),
+     pnorm(x,3,s=0) == c(0,0,0,0,  1 , 1,1,1,1,1),
+     dnorm(x,3,s=Inf) == 0,
+     pnorm(x,3,s=Inf) == c(0, rep(0.5, 8), 1))
> 
> ## 3 Test data from Wichura (1988) :
> stopifnot(
+     all.equal(qnorm(c( 0.25,  .001, 1e-20)),
+ 	      c(-0.6744897501960817, -3.090232306167814, -9.262340089798408),
+ 	      tolerance = 1e-15)
+   , ## extreme tail -- available on log scale only:
+     all.equal(qnorm(-1e5, log = TRUE), -447.1974945)
+ )
> 
> z <- rnorm(1000); all.equal(pnorm(z),  1 - pnorm(-z), tolerance = 1e-15)
[1] TRUE
> z <- c(-Inf,Inf,NA,NaN, rt(1000, df=2))
> z.ok <- z > -37.5 | !is.finite(z)
> for(df in 1:10) stopifnot(all.equal(pt(z, df), 1 - pt(-z,df), tolerance = 1e-15))
> 
> stopifnot(All.eq(pz <- pnorm(z), 1 - pnorm(z, lower=FALSE)),
+           All.eq(pz,		    pnorm(-z, lower=FALSE)),
+           All.eq(log(pz[z.ok]), pnorm(z[z.ok], log=TRUE)))
> y <- seq(-70,0, by = 10)
> cbind(y, "log(pnorm(y))"= log(pnorm(y)), "pnorm(y, log=T)"= pnorm(y, log=TRUE))
       y log(pnorm(y)) pnorm(y, log=T)
[1,] -70          -Inf   -2455.1676378
[2,] -60          -Inf   -1805.0135607
[3,] -50          -Inf   -1254.8313611
[4,] -40          -Inf    -804.6084420
[5,] -30  -454.3212440    -454.3212440
[6,] -20  -203.9171554    -203.9171554
[7,] -10   -53.2312852     -53.2312852
[8,]   0    -0.6931472      -0.6931472
> y <- c(1:15, seq(20,40, by=5))
> cbind(y, "log(pnorm(y))"= log(pnorm(y)), "pnorm(y, log=T)"= pnorm(y, log=TRUE),
+       "log(pnorm(-y))"= log(pnorm(-y)), "pnorm(-y, log=T)"= pnorm(-y, log=TRUE))
       y log(pnorm(y)) pnorm(y, log=T) log(pnorm(-y)) pnorm(-y, log=T)
 [1,]  1 -1.727538e-01   -1.727538e-01      -1.841022        -1.841022
 [2,]  2 -2.301291e-02   -2.301291e-02      -3.783184        -3.783184
 [3,]  3 -1.350810e-03   -1.350810e-03      -6.607726        -6.607726
 [4,]  4 -3.167174e-05   -3.167174e-05     -10.360101       -10.360101
 [5,]  5 -2.866516e-07   -2.866516e-07     -15.064998       -15.064998
 [6,]  6 -9.865877e-10   -9.865876e-10     -20.736769       -20.736769
 [7,]  7 -1.279865e-12   -1.279813e-12     -27.384307       -27.384307
 [8,]  8 -6.661338e-16   -6.220961e-16     -35.013437       -35.013437
 [9,]  9  0.000000e+00   -1.128588e-19     -43.628149       -43.628149
[10,] 10  0.000000e+00   -7.619853e-24     -53.231285       -53.231285
[11,] 11  0.000000e+00   -1.910660e-28     -63.824934       -63.824934
[12,] 12  0.000000e+00   -1.776482e-33     -75.410673       -75.410673
[13,] 13  0.000000e+00   -6.117164e-39     -87.989720       -87.989720
[14,] 14  0.000000e+00   -7.793537e-45    -101.563034      -101.563034
[15,] 15  0.000000e+00   -3.670966e-51    -116.131385      -116.131385
[16,] 20  0.000000e+00   -2.753624e-89    -203.917155      -203.917155
[17,] 25  0.000000e+00  -3.056697e-138    -316.639408      -316.639408
[18,] 30  0.000000e+00  -4.906714e-198    -454.321244      -454.321244
[19,] 35  0.000000e+00  -1.124911e-268    -616.975101      -616.975101
[20,] 40  0.000000e+00    0.000000e+00           -Inf      -804.608442
> ## Symmetry:
> y <- c(1:50,10^c(3:10,20,50,150,250))
> y <- c(-y,0,y)
> for(L in c(FALSE,TRUE))
+     stopifnot(identical(pnorm(-y, log= L),
+ 			pnorm(+y, log= L, lower=FALSE)))
> 
> ## Log norm
> stopifnot(All.eq(pz, plnorm(exp(z))))
> 
> 
> ###==========  p <-> q	Inversion consistency =====================
> ok <- 1e-5 < pz & pz < 1 - 1e-5
> all.equal(z[ok], qnorm(pz[ok]), tolerance = 1e-12)
[1] TRUE
> 
> ###===== Random numbers -- first, just output:
> 
> set.seed(123)
> # .Random.seed <- c(0L, 17292L, 29447L, 24113L)
> n <- 20
> ## for(pre in PDQR) { n <- paste("r",pre,sep=""); cat(n,": "); str(get(n))}
> (Rbeta	  <- rbeta    (n, shape1 = .8, shape2 = 2) )
 [1] 0.554206761 0.387924405 0.012541339 0.257889595 0.236064413 0.008248085
 [7] 0.136510774 0.618501837 0.028641308 0.151670292 0.242835922 0.551824427
[13] 0.127798688 0.087335901 0.509096247 0.382121566 0.773494885 0.639404676
[19] 0.066559813 0.227487378
> (Rbinom	  <- rbinom   (n, size = 55, prob = pi/16) )
 [1]  7 13 15 10 12  7 10  9 13 10 13 13 13 10 13 12 12  3 11  8
> (Rcauchy  <- rcauchy  (n, location = 12, scale = 2) )
 [1] 17.042930  6.592944 15.980645 12.728113 13.921385  8.570544 19.557102
 [8] 10.430404 12.669798 21.641273 11.905411 11.301537 11.254793 13.226057
[15] 12.870672  8.167569 15.735143  8.272588 15.159679 13.338095
> (Rchisq	  <- rchisq   (n, df = 3) )
 [1] 3.8641030 1.8367371 1.2291085 1.9151780 2.5901414 0.4522238 0.1018120
 [8] 0.4323865 2.4551985 2.8313420 1.1215175 8.6109152 0.2215425 2.6531221
[15] 2.9968656 2.1074132 4.2698236 8.2015930 1.1384935 0.6041294
> (Rexp	  <- rexp     (n, rate = 2) )
 [1] 0.719726313 0.232135206 0.134439025 0.731650280 0.768696989 0.002299563
 [7] 0.554382734 0.149985159 0.596001504 0.557464352 0.033687945 0.240334360
[13] 0.785227170 0.129973053 0.928461144 0.231609810 0.118017873 0.591049711
[19] 0.029835686 0.201619221
> (Rf	  <- rf	      (n, df1 = 12, df2 = 6) )
 [1] 0.7630651 0.4158576 0.4052799 0.8862771 2.0435073 2.5383554 0.8259689
 [8] 1.5831584 0.6362728 0.7221707 1.1665204 0.4821276 1.1985013 0.8955592
[15] 1.5186453 0.7333794 0.3795005 0.7464008 1.5543317 5.8525268
> (Rgamma	  <- rgamma   (n, shape = 2, scale = 5) )
 [1] 11.843056  6.993596 12.413160  5.980309  5.970495 15.674812  7.382163
 [8] 11.343789  5.392676  3.766812 15.331767  2.451095 13.046780 10.999054
[15]  9.929606  6.239055  3.588397 26.196937  3.647236 23.051571
> (Rgeom	  <- rgeom    (n, prob = pi/16) )
 [1] 8 2 4 2 3 5 5 4 1 8 1 6 1 6 8 3 8 7 5 0
> (Rhyper	  <- rhyper   (n, m = 40, n = 30, k = 20) )
 [1] 12 11 12 10 13 11 13 13  8 12 12 15 11 11 11 12  9 11 12 12
> (Rlnorm	  <- rlnorm   (n, meanlog = -1, sdlog = 3) )
 [1] 1.873201e-02 5.609747e+01 9.793124e-02 4.203730e-02 9.015415e-03
 [6] 7.795979e-03 6.574837e-02 2.348924e+00 1.027326e+01 3.073340e+00
[11] 1.235670e-01 4.401014e-01 4.443226e-02 4.278129e-02 5.227596e+00
[16] 1.747860e-02 1.297853e+02 2.805625e-01 7.002040e-01 4.013193e-02
> (Rlogis	  <- rlogis   (n, location = 12, scale = 2) )
 [1] 10.139287 14.905198  7.466552 15.060325 11.415297 14.260524 13.347886
 [8] 11.554385 13.039986 -3.344515  9.436395 13.741216  9.411881 14.951552
[15] 10.378811 13.579174 17.257244  7.933965  8.157203 13.491243
> (Rnbinom  <- rnbinom  (n, size = 7, prob = .01) )
 [1]  563  315  519  707  614  637  560 1086 1088  842  522  787  576  673  607
[16]  443 1023  590  663  947
> (Rnorm	  <- rnorm    (n, mean = -1, sd = 3) )
 [1] -0.3732778 -5.6447208 -4.5214067  2.5364456 -2.2579283 -6.5577259
 [7]  0.4577794 -2.4157378  1.5402158  1.7773118 -1.9624917 -4.9029209
[13]  6.3787040  1.3350950  1.0696427 -6.1133966  2.3161164 -1.9415584
[19] -2.0400782  2.1153214
> (Rpois	  <- rpois    (n, lambda = 12) )
 [1] 10 15 14  9 15 12 16  7  7 14 11 14 16 12 15  7  7 13 15 12
> (Rsignrank<- rsignrank(n, n = 47) )
 [1] 531 522 385 583 731 461 542 590 485 582 571 475 702 537 395 571 603 515 477
[20] 609
> (Rt	  <- rt	      (n, df = 11) )
 [1] -0.499846093 -1.670698483 -0.729757704  0.009882459  1.518830804
 [6] -1.105480914 -0.073913420 -0.323601935 -0.309210484  1.071490670
[11] -1.730648357 -0.283141485  0.205309956  0.409676847  3.137494386
[16]  0.187070891 -1.154130926  0.970408359 -1.340678412  1.437769232
> ## Rt2 below (to preserve the following random numbers!)
> (Runif	  <- runif    (n, min = .2, max = 2) )
 [1] 1.3077353 0.6915398 0.3146723 0.3301459 1.7723202 0.8530850 1.9021199
 [8] 1.2049741 1.1227901 1.9059890 0.3976429 1.6740336 0.8792158 1.2019539
[15] 0.7649681 0.3482172 1.6038876 0.5478225 1.3800538 0.9832647
> (Rweibull <- rweibull (n, shape = 3, scale = 2) )
 [1] 2.3507552 0.3498804 1.8610901 3.0250054 0.8043691 1.0338586 1.5966081
 [8] 2.1295946 1.1916615 0.7320637 1.0860054 2.5452588 0.6323843 1.6726976
[15] 3.4711560 2.5185134 3.4519830 1.1740161 1.2526410 0.8298120
> (Rwilcox  <- rwilcox  (n, m = 13, n = 17) )
 [1]  80  92 104  92 116  55 115 105 102 133  98 116 149 118  93  73 126 126 124
[20]  99
> (Rt2	  <- rt	      (n, df = 1.01))
 [1] -0.4507473  7.9567158  1.9317506 -4.5681636 -0.7939303 -1.3143350
 [7] -0.2751314 -1.2305970 -1.2321505  2.2912240  4.3743122 -1.3400950
[13] -0.4901184  0.7864050  0.9570928 50.6881882  0.2577983  0.2478237
[19]  1.3330087  0.3807687
> 
> (Pbeta	  <- pbeta    (Rbeta, shape1 = .8, shape2 = 2) )
 [1] 0.84605840 0.69836904 0.05389148 0.53895245 0.50764485 0.03861536
 [7] 0.34373485 0.88868811 0.10358645 0.37126511 0.51751532 0.84433212
[13] 0.32741433 0.24605563 0.81153016 0.69214938 0.96181357 0.90094249
[19] 0.19989455 0.49493375
> (Pbinom	  <- pbinom   (Rbinom, size = 55, prob = pi/16) )
 [1] 0.128811009 0.822082686 0.939887897 0.473358125 0.725714645 0.128811009
 [7] 0.473358125 0.340181320 0.822082686 0.473358125 0.822082686 0.822082686
[13] 0.822082686 0.473358125 0.822082686 0.725714645 0.725714645 0.002917878
[19] 0.606468256 0.221684127
> (Pcauchy  <- pcauchy  (Rcauchy, location = 12, scale = 2) )
 [1] 0.8798165 0.1127710 0.8517979 0.6111354 0.7436195 0.1680556 0.9176468
 [8] 0.2881958 0.6028646 0.9348927 0.4849570 0.3930511 0.3864691 0.6750527
[15] 0.6306957 0.1531019 0.8435165 0.1567581 0.8203732 0.6876911
> (Pchisq	  <- pchisq   (Rchisq, df = 3) )
 [1] 0.723482597 0.393026646 0.253968309 0.409802926 0.540779498 0.070744975
 [7] 0.008380908 0.066526148 0.516557851 0.581633585 0.228117628 0.965062738
[13] 0.025961119 0.551747428 0.607891302 0.449585516 0.766236467 0.957975966
[19] 0.232209338 0.104513603
> (Pexp	  <- pexp     (Rexp, rate = 2) )
 [1] 0.762942518 0.371406447 0.235763552 0.768528971 0.785059488 0.004588567
 [7] 0.670033937 0.259159790 0.696387483 0.672061341 0.065156264 0.381630263
[13] 0.792049319 0.228906859 0.843847516 0.370745580 0.210247557 0.693365692
[19] 0.057925925 0.331847235
> (Pf	  <- pf	      (Rf, df1 = 12, df2 = 6) )
 [1] 0.32410218 0.09232470 0.08630331 0.40259602 0.80440198 0.86907389
 [7] 0.36505563 0.70322821 0.23758524 0.29664970 0.55170692 0.13291359
[13] 0.56607937 0.40821091 0.68428769 0.30422992 0.07227779 0.31298492
[19] 0.69494050 0.97980295
> (Pgamma	  <- pgamma   (Rgamma, shape = 2, scale = 5) )
 [1] 0.6846621 0.4077251 0.7091198 0.3359489 0.3352389 0.8201234 0.4342559
 [8] 0.6618778 0.2931040 0.1745504 0.8105456 0.0872524 0.7344167 0.6453838
[15] 0.5901700 0.3545801 0.1619732 0.9669094 0.1661030 0.9441847
> (Pgeom	  <- pgeom    (Rgeom, prob = pi/16) )
 [1] 0.8601686 0.4809591 0.6647753 0.4809591 0.5828725 0.7305965 0.7305965
 [8] 0.6647753 0.3541459 0.8601686 0.3541459 0.7834938 0.3541459 0.7834938
[15] 0.8601686 0.5828725 0.8601686 0.8260047 0.7305965 0.1963495
> (Phyper	  <- phyper   (Rhyper, m = 40, n = 30, k = 20) )
 [1] 0.71494883 0.51295659 0.71494883 0.30864260 0.86627413 0.51295659
 [7] 0.86627413 0.86627413 0.05904998 0.71494883 0.71494883 0.98680472
[13] 0.51295659 0.51295659 0.51295659 0.71494883 0.15132082 0.51295659
[19] 0.71494883 0.71494883
> (Plnorm	  <- plnorm   (Rlnorm, meanlog = -1, sdlog = 3) )
 [1] 0.16047510 0.95310121 0.32954742 0.23481968 0.10817853 0.09944582
 [7] 0.28299287 0.73170766 0.86646776 0.76039954 0.35805697 0.52382260
[13] 0.24053070 0.23661975 0.81182742 0.15491172 0.97472581 0.46401663
[19] 0.58493655 0.23009692
> (Plogis	  <- plogis   (Rlogis, location = 12, scale = 2) )
 [1] 0.2828523996 0.8103980662 0.0939166071 0.8220300712 0.4274282807
 [6] 0.7558872597 0.6623855066 0.4445273969 0.6271461844 0.0004653491
[11] 0.2172435140 0.7048722457 0.2151664882 0.8139337213 0.3077638743
[16] 0.6877426857 0.9326810790 0.1157796648 0.1277056879 0.6782238574
> (Pnbinom  <- pnbinom  (Rnbinom, size = 7, prob = .01) )
 [1] 0.34539127 0.04498289 0.27616523 0.57177025 0.42755552 0.46435569
 [7] 0.34059007 0.91950991 0.92034921 0.74409223 0.28078034 0.68045166
[13] 0.36627489 0.52069882 0.41628414 0.16798077 0.88887110 0.38885301
[19] 0.50524696 0.83942274
> (Pnorm	  <- pnorm    (Rnorm, mean = -1, sd = 3) )
 [1] 0.58273974 0.06078223 0.12023713 0.88076412 0.33749500 0.03197163
 [7] 0.68649035 0.31849459 0.80142958 0.82271739 0.37416939 0.09663374
[13] 0.99304478 0.78182307 0.75486546 0.04414707 0.86550081 0.37681641
[19] 0.36441108 0.85046748
> (Ppois	  <- ppois    (Rpois, lambda = 12) )
 [1] 0.3472294 0.8444157 0.7720245 0.2423922 0.8444157 0.5759652 0.8987090
 [8] 0.0895045 0.0895045 0.7720245 0.4615973 0.7720245 0.8987090 0.5759652
[15] 0.8444157 0.0895045 0.0895045 0.6815356 0.8444157 0.5759652
> (Psignrank<- psignrank(Rsignrank, n = 47) )
 [1] 0.36663746 0.33169387 0.02921581 0.58098562 0.96191228 0.14044836
 [7] 0.41083885 0.60942723 0.20480997 0.57688473 0.53132960 0.17614726
[13] 0.92806867 0.39057277 0.03720129 0.53132960 0.66065154 0.30545809
[19] 0.18167093 0.68340768
> (Pt	  <- pt	      (Rt,  df = 11) )
 [1] 0.31351556 0.06147847 0.24039508 0.50385399 0.92149468 0.14627181
 [7] 0.47120313 0.37615421 0.38146974 0.84655256 0.05571530 0.39116296
[13] 0.57946003 0.65504612 0.99527474 0.57249461 0.13645270 0.82365134
[19] 0.10352754 0.91083331
> (Pt2	  <- pt	      (Rt2, df = 1.01) )
 [1] 0.36490878 0.96094812 0.84890254 0.06767383 0.28587141 0.20622746
 [7] 0.41436382 0.21644090 0.21624373 0.86999347 0.92939659 0.20324856
[13] 0.35461982 0.71264902 0.74365588 0.99395091 0.58047323 0.57748238
[19] 0.79593933 0.61605276
> (Punif	  <- punif    (Runif, min = .2, max = 2) )
 [1] 0.61540849 0.27307767 0.06370684 0.07230328 0.87351121 0.36282500
 [7] 0.94562215 0.55831894 0.51266115 0.94777166 0.10980160 0.81890755
[13] 0.37734211 0.55664105 0.31387119 0.08234288 0.77993754 0.19323473
[19] 0.65558547 0.43514705
> (Pweibull <- pweibull (Rweibull, shape = 3, scale = 2) )
 [1] 0.802851673 0.005339578 0.553257134 0.968573336 0.062983442 0.129016046
 [7] 0.398753717 0.700984544 0.190653725 0.047857615 0.147946029 0.872690264
[13] 0.031117650 0.442898991 0.994635553 0.864236080 0.994152898 0.183125938
[19] 0.217836412 0.068933714
> (Pwilcox  <- pwilcox  (Rwilcox, m = 13, n = 17) )
 [1] 0.106381913 0.228522998 0.402346437 0.228522998 0.597653563 0.009833454
 [7] 0.581622297 0.418377703 0.370781167 0.829464624 0.310291237 0.597653563
[13] 0.948370460 0.629218833 0.241282199 0.061358235 0.745593076 0.745593076
[19] 0.718298294 0.325029106
> 
> dbeta	 (Rbeta, shape1 = .8, shape2 = 2)
 [1] 0.7223732 1.0651681 3.4136000 1.4013424 1.4682944 3.7281631 1.8517657
 [8] 0.6047651 2.8467121 1.7813410 1.4470712 0.7268595 1.8952825 2.1401041
[15] 0.8090938 1.0785126 0.3433595 0.5678419 2.3110352 1.4958106
> dbinom	 (Rbinom, size = 55, prob = pi/16)
 [1] 0.063354298 0.096368041 0.047170694 0.133176805 0.119246389 0.063354298
 [7] 0.133176805 0.118497194 0.096368041 0.133176805 0.096368041 0.096368041
[13] 0.096368041 0.133176805 0.096368041 0.119246389 0.119246389 0.002298608
[19] 0.133110131 0.092873118
> dcauchy	 (Rcauchy, location = 12, scale = 2)
 [1] 0.021630824 0.019154376 0.032078738 0.140529537 0.082766881 0.040391664
 [7] 0.010417640 0.098492589 0.143104718 0.006566194 0.158799748 0.141854068
[13] 0.139752618 0.115681430 0.133797911 0.034066567 0.035463730 0.035578069
[19] 0.045526265 0.109942144
> dchisq	 (Rchisq, df = 3)
 [1] 0.11359392 0.21581974 0.23922560 0.21190393 0.17584488 0.21398711
 [7] 0.12097660 0.21132681 0.18315290 0.16296267 0.24114522 0.01579786
[13] 0.16808590 0.17245292 0.15434146 0.20191378 0.09748429 0.01891934
[19] 0.24090989 0.22923980
> dexp	 (Rexp, rate = 2)
 [1] 0.4741150 1.2571871 1.5284729 0.4629421 0.4298810 1.9908229 0.6599321
 [8] 1.4816804 0.6072250 0.6558773 1.8696875 1.2367395 0.4159014 1.5421863
[15] 0.3123050 1.2585088 1.5795049 0.6132686 1.8841481 1.3363055
> df	 (Rf, df1 = 12, df2 = 6)
 [1] 0.664042111 0.576090730 0.562274916 0.607318497 0.167959027 0.100284610
 [7] 0.637144174 0.282709564 0.693573655 0.678003283 0.457493848 0.643454728
[13] 0.441386316 0.602499118 0.304762080 0.674516557 0.524996364 0.670131912
[19] 0.292349736 0.008558897
> dgamma	 (Rgamma, shape = 2, scale = 5)
 [1] 0.044345445 0.069072391 0.041471424 0.072333682 0.072356861 0.027275042
 [7] 0.067458566 0.046936330 0.073360525 0.070933729 0.028572738 0.060051063
[13] 0.038400284 0.048758423 0.054515175 0.071657454 0.070028735 0.005557423
[19] 0.070344309 0.009173320
> dgeom	 (Rgeom, prob = pi/16)
 [1] 0.03416390 0.12681315 0.08190279 0.12681315 0.10191344 0.06582121
 [7] 0.06582121 0.08190279 0.15779640 0.03416390 0.15779640 0.05289725
[13] 0.15779640 0.05289725 0.03416390 0.10191344 0.03416390 0.04251090
[19] 0.06582121 0.19634954
> dhyper	 (Rhyper, m = 40, n = 30, k = 20)
 [1] 0.20199224 0.20431399 0.20199224 0.15732178 0.15132529 0.20431399
 [7] 0.15132529 0.15132529 0.04108936 0.20199224 0.20199224 0.03541012
[13] 0.20431399 0.20431399 0.20431399 0.20199224 0.09227084 0.20431399
[19] 0.20199224 0.20199224
> dlnorm	 (Rlnorm, meanlog = -1, sdlog = 3)
 [1] 4.3380954158 0.0005822436 1.2319845653 2.4357017049 6.8694422907
 [6] 7.4733929149 1.7154009417 0.0467724227 0.0069920836 0.0336865537
[11] 1.0073244481 0.3016205039 2.3349967070 2.4034440204 0.0172006936
[16] 4.5426601250 0.0001514856 0.4720497698 0.1855964185 2.5226771709
> dlogis	 (Rlogis, location = 12, scale = 2)
 [1] 0.1014234598 0.0768265203 0.0425481390 0.0731483166 0.1223666728
 [6] 0.0922608552 0.1118154736 0.1234613952 0.1169169239 0.0002325663
[11] 0.0850243848 0.1040136815 0.0844349353 0.0757228093 0.1065226360
[16] 0.1073763420 0.0313935419 0.0511873670 0.0556984726 0.1091181283
> dnbinom	 (Rnbinom, size = 7, prob = .01)
 [1] 0.0016012889 0.0006114669 0.0015340423 0.0014659766 0.0016087509
 [6] 0.0015899813 0.0015985567 0.0004225217 0.0004186957 0.0010719579
[11] 0.0015404998 0.0012444310 0.0016100854 0.0015372764 0.0016118864
[16] 0.0012822374 0.0005567006 0.0016141573 0.0015545113 0.0007531913
> dnorm	 (Rnorm, mean = -1, sd = 3)
 [1] 0.130110398 0.040112199 0.066772988 0.066380399 0.121789513 0.023907371
 [7] 0.118172118 0.118967900 0.092918367 0.086632832 0.126309889 0.057050581
[13] 0.006458952 0.098226734 0.104819164 0.031112294 0.072188495 0.126589892
[19] 0.125224300 0.077557930
> dpois	 (Rpois, lambda = 12)
 [1] 0.10483726 0.07239112 0.09048890 0.08736438 0.07239112 0.11436792
 [7] 0.05429334 0.04368219 0.04368219 0.09048890 0.11436792 0.09048890
[13] 0.05429334 0.11436792 0.07239112 0.04368219 0.04368219 0.10557038
[19] 0.07239112 0.11436792
> dsignrank(Rsignrank, n = 47)
 [1] 0.0039429676 0.0038018425 0.0007151971 0.0041008982 0.0009032331
 [6] 0.0023508267 0.0040737886 0.0040318201 0.0029828933 0.0041090875
[11] 0.0041704554 0.0027222983 0.0014781510 0.0040203094 0.0008698303
[16] 0.0041704554 0.0038520708 0.0036731856 0.0027749926 0.0037486048
> dt	 (Rt, df = 11)
 [1] 0.340823726 0.100413165 0.293668976 0.389968983 0.124439520 0.207270198
 [7] 0.388829635 0.368437694 0.370255866 0.214961569 0.091948602 0.373362745
[13] 0.381142118 0.356119169 0.008424401 0.382627645 0.196428441 0.238239645
[19] 0.157280563 0.138777161
> dunif	 (Runif, min = .2, max = 2)
 [1] 0.5555556 0.5555556 0.5555556 0.5555556 0.5555556 0.5555556 0.5555556
 [8] 0.5555556 0.5555556 0.5555556 0.5555556 0.5555556 0.5555556 0.5555556
[15] 0.5555556 0.5555556 0.5555556 0.5555556 0.5555556 0.5555556
> dweibull (Rweibull, shape = 3, scale = 2)
 [1] 0.40854433 0.04566100 0.58026142 0.10784049 0.22734704 0.34911114
 [7] 0.57475176 0.50853257 0.43099423 0.19135106 0.37684466 0.30928353
[13] 0.14529960 0.58452097 0.02423843 0.32292685 0.02612818 0.42221582
[19] 0.46023762 0.24042036
> dwilcox	 (Rwilcox, m = 13, n = 17)
 [1] 0.007553976 0.012383032 0.015873676 0.012383032 0.016031266 0.001126580
 [7] 0.016163798 0.016031266 0.015486075 0.010799421 0.014448507 0.016031266
[13] 0.004704515 0.015691594 0.012759201 0.005024180 0.013477697 0.013477697
[19] 0.014141025 0.014737869
> 
> ## Check q*(p*(.)) = identity
> All.eq(Rbeta,	  qbeta	   (Pbeta, shape1 = .8, shape2 = 2))
[1] TRUE
> All.eq(Rbinom,	  qbinom   (Pbinom, size = 55, prob = pi/16))
[1] TRUE
> All.eq(Rcauchy,	  qcauchy  (Pcauchy, location = 12, scale = 2))
[1] TRUE
> All.eq(Rchisq,	  qchisq   (Pchisq, df = 3))
[1] TRUE
> All.eq(Rexp,	  qexp	   (Pexp, rate = 2))
[1] TRUE
> All.eq(Rf,	  qf	   (Pf, df1 = 12, df2 = 6))
[1] TRUE
> All.eq(Rgamma,	  qgamma   (Pgamma, shape = 2, scale = 5))
[1] TRUE
> All.eq(Rgeom,	  qgeom	   (Pgeom, prob = pi/16))
[1] TRUE
> All.eq(Rhyper,	  qhyper   (Phyper, m = 40, n = 30, k = 20))
[1] TRUE
> All.eq(Rlnorm,	  qlnorm   (Plnorm, meanlog = -1, sdlog = 3))
[1] TRUE
> All.eq(Rlogis,	  qlogis   (Plogis, location = 12, scale = 2))
[1] TRUE
> All.eq(Rnbinom,	  qnbinom  (Pnbinom, size = 7, prob = .01))
[1] TRUE
> All.eq(Rnorm,	  qnorm	   (Pnorm, mean = -1, sd = 3))
[1] TRUE
> All.eq(Rpois,	  qpois	   (Ppois, lambda = 12))
[1] TRUE
> All.eq(Rsignrank, qsignrank(Psignrank, n = 47))
[1] TRUE
> All.eq(Rt,	  qt	   (Pt,	 df = 11))
[1] TRUE
> All.eq(Rt2,	  qt	   (Pt2, df = 1.01))
[1] TRUE
> All.eq(Runif,	  qunif	   (Punif, min = .2, max = 2))
[1] TRUE
> All.eq(Rweibull,  qweibull (Pweibull, shape = 3, scale = 2))
[1] TRUE
> All.eq(Rwilcox,	  qwilcox  (Pwilcox, m = 13, n = 17))
[1] TRUE
> 
> ## Same with "upper tail":
> All.eq(Rbeta,	  qbeta	   (1- Pbeta, shape1 = .8, shape2 = 2, lower=F))
[1] TRUE
> All.eq(Rbinom,	  qbinom   (1- Pbinom, size = 55, prob = pi/16, lower=F))
[1] TRUE
> All.eq(Rcauchy,	  qcauchy  (1- Pcauchy, location = 12, scale = 2, lower=F))
[1] TRUE
> All.eq(Rchisq,	  qchisq   (1- Pchisq, df = 3, lower=F))
[1] TRUE
> All.eq(Rexp,	  qexp	   (1- Pexp, rate = 2, lower=F))
[1] TRUE
> All.eq(Rf,	  qf	   (1- Pf, df1 = 12, df2 = 6, lower=F))
[1] TRUE
> All.eq(Rgamma,	  qgamma   (1- Pgamma, shape = 2, scale = 5, lower=F))
[1] TRUE
> All.eq(Rgeom,	  qgeom	   (1- Pgeom, prob = pi/16, lower=F))
[1] TRUE
> All.eq(Rhyper,	  qhyper   (1- Phyper, m = 40, n = 30, k = 20, lower=F))
[1] TRUE
> All.eq(Rlnorm,	  qlnorm   (1- Plnorm, meanlog = -1, sdlog = 3, lower=F))
[1] TRUE
> All.eq(Rlogis,	  qlogis   (1- Plogis, location = 12, scale = 2, lower=F))
[1] TRUE
> All.eq(Rnbinom,	  qnbinom  (1- Pnbinom, size = 7, prob = .01, lower=F))
[1] TRUE
> All.eq(Rnorm,	  qnorm	   (1- Pnorm, mean = -1, sd = 3,lower=F))
[1] TRUE
> All.eq(Rpois,	  qpois	   (1- Ppois, lambda = 12, lower=F))
[1] TRUE
> All.eq(Rsignrank, qsignrank(1- Psignrank, n = 47, lower=F))
[1] TRUE
> All.eq(Rt,	  qt	   (1- Pt,  df = 11,   lower=F))
[1] TRUE
> All.eq(Rt2,	  qt	   (1- Pt2, df = 1.01, lower=F))
[1] TRUE
> All.eq(Runif,	  qunif	   (1- Punif, min = .2, max = 2, lower=F))
[1] TRUE
> All.eq(Rweibull,  qweibull (1- Pweibull, shape = 3, scale = 2, lower=F))
[1] TRUE
> All.eq(Rwilcox,	  qwilcox  (1- Pwilcox, m = 13, n = 17, lower=F))
[1] TRUE
> 
> ## Check q*(p* ( log ), log) = identity
> All.eq(Rbeta,	  qbeta	   (log(Pbeta), shape1 = .8, shape2 = 2, log=TRUE))
[1] TRUE
> All.eq(Rbinom,	  qbinom   (log(Pbinom), size = 55, prob = pi/16, log=TRUE))
[1] TRUE
> All.eq(Rcauchy,	  qcauchy  (log(Pcauchy), location = 12, scale = 2, log=TRUE))
[1] TRUE
> All.eq(Rchisq,    qchisq   (log(Pchisq), df = 3, log=TRUE))
[1] TRUE
> All.eq(Rexp,	  qexp	   (log(Pexp), rate = 2, log=TRUE))
[1] TRUE
> All.eq(Rf,	  qf	   (log(Pf), df1= 12, df2= 6, log=TRUE))
[1] TRUE
> All.eq(Rgamma,	  qgamma   (log(Pgamma), shape = 2, scale = 5, log=TRUE))
[1] TRUE
> All.eq(Rgeom,	  qgeom	   (log(Pgeom), prob = pi/16, log=TRUE))
[1] TRUE
> All.eq(Rhyper,	  qhyper   (log(Phyper), m = 40, n = 30, k = 20, log=TRUE))
[1] TRUE
> All.eq(Rlnorm,	  qlnorm   (log(Plnorm), meanlog = -1, sdlog = 3, log=TRUE))
[1] TRUE
> All.eq(Rlogis,	  qlogis   (log(Plogis), location = 12, scale = 2, log=TRUE))
[1] TRUE
> All.eq(Rnbinom,	  qnbinom  (log(Pnbinom), size = 7, prob = .01, log=TRUE))
[1] TRUE
> All.eq(Rnorm,	  qnorm	   (log(Pnorm), mean = -1, sd = 3, log=TRUE))
[1] TRUE
> All.eq(Rpois,	  qpois	   (log(Ppois), lambda = 12, log=TRUE))
[1] TRUE
> All.eq(Rsignrank, qsignrank(log(Psignrank), n = 47, log=TRUE))
[1] TRUE
> All.eq(Rt,	  qt	   (log(Pt), df = 11, log=TRUE))
[1] TRUE
> All.eq(Rt2,	  qt	   (log(Pt2), df = 1.01, log=TRUE))
[1] TRUE
> All.eq(Runif,	  qunif	   (log(Punif), min = .2, max = 2, log=TRUE))
[1] TRUE
> All.eq(Rweibull,  qweibull (log(Pweibull), shape = 3, scale = 2, log=TRUE))
[1] TRUE
> All.eq(Rwilcox,	  qwilcox  (log(Pwilcox), m = 13, n = 17, log=TRUE))
[1] TRUE
> 
> ## same q*(p* (log) log) with upper tail:
> 
> All.eq(Rbeta,	  qbeta	   (log1p(-Pbeta), shape1 = .8, shape2 = 2, lower=F, log=T))
[1] TRUE
> All.eq(Rbinom,	  qbinom   (log1p(-Pbinom), size = 55, prob = pi/16, lower=F, log=T))
[1] TRUE
> All.eq(Rcauchy,	  qcauchy  (log1p(-Pcauchy), location = 12, scale = 2, lower=F, log=T))
[1] TRUE
> All.eq(Rchisq,	  qchisq   (log1p(-Pchisq), df = 3, lower=F, log=T))
[1] TRUE
> All.eq(Rexp,	  qexp	   (log1p(-Pexp), rate = 2, lower=F, log=T))
[1] TRUE
> All.eq(Rf,	  qf	   (log1p(-Pf), df1 = 12, df2 = 6, lower=F, log=T))
[1] TRUE
> All.eq(Rgamma,	  qgamma   (log1p(-Pgamma), shape = 2, scale = 5, lower=F, log=T))
[1] TRUE
> All.eq(Rgeom,	  qgeom	   (log1p(-Pgeom), prob = pi/16, lower=F, log=T))
[1] TRUE
> All.eq(Rhyper,	  qhyper   (log1p(-Phyper), m = 40, n = 30, k = 20, lower=F, log=T))
[1] TRUE
> All.eq(Rlnorm,	  qlnorm   (log1p(-Plnorm), meanlog = -1, sdlog = 3, lower=F, log=T))
[1] TRUE
> All.eq(Rlogis,	  qlogis   (log1p(-Plogis), location = 12, scale = 2, lower=F, log=T))
[1] TRUE
> All.eq(Rnbinom,	  qnbinom  (log1p(-Pnbinom), size = 7, prob = .01, lower=F, log=T))
[1] TRUE
> All.eq(Rnorm,	  qnorm	   (log1p(-Pnorm), mean = -1, sd = 3, lower=F, log=T))
[1] TRUE
> All.eq(Rpois,	  qpois	   (log1p(-Ppois), lambda = 12, lower=F, log=T))
[1] TRUE
> All.eq(Rsignrank, qsignrank(log1p(-Psignrank), n = 47, lower=F, log=T))
[1] TRUE
> All.eq(Rt,	  qt	   (log1p(-Pt ), df = 11,   lower=F, log=T))
[1] TRUE
> All.eq(Rt2,	  qt	   (log1p(-Pt2), df = 1.01, lower=F, log=T))
[1] TRUE
> All.eq(Runif,	  qunif	   (log1p(-Punif), min = .2, max = 2, lower=F, log=T))
[1] TRUE
> All.eq(Rweibull,  qweibull (log1p(-Pweibull), shape = 3, scale = 2, lower=F, log=T))
[1] TRUE
> All.eq(Rwilcox,	  qwilcox  (log1p(-Pwilcox), m = 13, n = 17, lower=F, log=T))
[1] TRUE
> 
> 
> ## Check log( upper.tail ):
> All.eq(log1p(-Pbeta),	  pbeta	   (Rbeta, shape1 = .8, shape2 = 2, lower=F, log=T))
[1] TRUE
> All.eq(log1p(-Pbinom),	  pbinom   (Rbinom, size = 55, prob = pi/16, lower=F, log=T))
[1] TRUE
> All.eq(log1p(-Pcauchy),	  pcauchy  (Rcauchy, location = 12, scale = 2, lower=F, log=T))
[1] TRUE
> All.eq(log1p(-Pchisq),	  pchisq   (Rchisq, df = 3, lower=F, log=T))
[1] TRUE
> All.eq(log1p(-Pexp),	  pexp	   (Rexp, rate = 2, lower=F, log=T))
[1] TRUE
> All.eq(log1p(-Pf),	  pf	   (Rf, df1 = 12, df2 = 6, lower=F, log=T))
[1] TRUE
> All.eq(log1p(-Pgamma),	  pgamma   (Rgamma, shape = 2, scale = 5, lower=F, log=T))
[1] TRUE
> All.eq(log1p(-Pgeom),	  pgeom	   (Rgeom, prob = pi/16, lower=F, log=T))
[1] TRUE
> All.eq(log1p(-Phyper),	  phyper   (Rhyper, m = 40, n = 30, k = 20, lower=F, log=T))
[1] TRUE
> All.eq(log1p(-Plnorm),	  plnorm   (Rlnorm, meanlog = -1, sdlog = 3, lower=F, log=T))
[1] TRUE
> All.eq(log1p(-Plogis),	  plogis   (Rlogis, location = 12, scale = 2, lower=F, log=T))
[1] TRUE
> All.eq(log1p(-Pnbinom),	  pnbinom  (Rnbinom, size = 7, prob = .01, lower=F, log=T))
[1] TRUE
> All.eq(log1p(-Pnorm),	  pnorm	   (Rnorm, mean = -1, sd = 3, lower=F, log=T))
[1] TRUE
> All.eq(log1p(-Ppois),	  ppois	   (Rpois, lambda = 12, lower=F, log=T))
[1] TRUE
> All.eq(log1p(-Psignrank), psignrank(Rsignrank, n = 47, lower=F, log=T))
[1] TRUE
> All.eq(log1p(-Pt),	  pt	   (Rt, df = 11,   lower=F, log=T))
[1] TRUE
> All.eq(log1p(-Pt2),	  pt	   (Rt2,df = 1.01, lower=F, log=T))
[1] TRUE
> All.eq(log1p(-Punif),	  punif	   (Runif, min = .2, max = 2, lower=F, log=T))
[1] TRUE
> All.eq(log1p(-Pweibull),  pweibull (Rweibull, shape = 3, scale = 2, lower=F, log=T))
[1] TRUE
> All.eq(log1p(-Pwilcox),	  pwilcox  (Rwilcox, m = 13, n = 17, lower=F, log=T))
[1] TRUE
> 
> 
> ### (Extreme) tail tests added more recently:
> All.eq(1, -1e-17/ pexp(qexp(-1e-17, log=TRUE),log=TRUE))
[1] TRUE
> abs(pgamma(30,100, lower=FALSE, log=TRUE) + 7.3384686328784e-24) < 1e-36
[1] TRUE
> All.eq(1, pcauchy(-1e20)           /  3.18309886183791e-21)
[1] TRUE
> All.eq(1, pcauchy(+1e15, log=TRUE) / -3.18309886183791e-16)## PR#6756
[1] TRUE
> x <- 10^(ex <- c(1,2,5*(1:5),50,100,200,300,Inf))
> for(a in x[ex > 10]) ## improve pt() : cbind(x,t= pt(-x, df=1), C=pcauchy(-x))
+     stopifnot(all.equal(pt(-a, df=1), pcauchy(-a), tolerance = 1e-15))
> ## for PR#7902:
> ex <- -c(rev(1/x), ex)
> All.eq(-x, qcauchy(pcauchy(-x)))
[1] TRUE
> All.eq(+x, qcauchy(pcauchy(+x, log=TRUE), log=TRUE))
[1] TRUE
> All.eq(1/x, pcauchy(qcauchy(1/x)))
[1] TRUE
> All.eq(ex,  pcauchy(qcauchy(ex, log=TRUE), log=TRUE))
[1] TRUE
> II <- c(-Inf,Inf)
> stopifnot(pcauchy(II) == 0:1, qcauchy(0:1) == II,
+           pcauchy(II, log=TRUE) == c(-Inf,0),
+           qcauchy(c(-Inf,0), log=TRUE) == II)
> ## PR#15521 :
> p <- 1 - 1/4096
> stopifnot(all.equal(qcauchy(p), 1303.7970381453319163, tolerance = 1e-14))
> 
> pr <- 1e-23 ## PR#6757
> stopifnot(all.equal(pr^ 12, pbinom(11, 12, prob= pr,lower=FALSE),
+                     tolerance = 1e-12, scale= 1e-270))
> ## pbinom(.) gave 0 in R 1.9.0
> pp <- 1e-17 ## PR#6792
> stopifnot(all.equal(2*pp, pgeom(1, pp), scale= 1e-20))
> ## pgeom(.) gave 0 in R 1.9.0
> 
> x <- 10^(100:295)
> sapply(c(1e-250, 1e-25, 0.9, 1.1, 101, 1e10, 1e100),
+        function(shape)
+        All.eq(-x, pgamma(x, shape=shape, lower=FALSE, log=TRUE)))
[1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE
> x <- 2^(-1022:-900)
> ## where all completely off in R 2.0.1
> all.equal(pgamma(x, 10, log = TRUE) - 10*log(x),
+           rep(-15.104412573076, length(x)), tolerance = 1e-12)# 3.984e-14 (i386)
[1] TRUE
> all.equal(pgamma(x, 0.1, log = TRUE) - 0.1*log(x),
+           rep(0.0498724412598364, length(x)), tolerance = 1e-13)# 7e-16 (i386)
[1] TRUE
> 
> All.eq(dpois(  10*1:2, 3e-308, log=TRUE),
+        c(-7096.08037610806, -14204.2875435307))
[1] TRUE
> All.eq(dpois(1e20, 1e-290, log=TRUE), -7.12801378828154e+22)
[1] TRUE
> ## all gave -Inf in R 2.0.1
> 
> 
> ## Inf df in pf etc.
> # apparently pf(df2=Inf) worked in 2.0.1 (undocumented) but df did not.
> x <- c(1/pi, 1, pi)
> oo <- options(digits = 8)
> df(x, 3, 1e6)
[1] 0.72553184 0.46254030 0.03300701
> df(x, 3, Inf)
[1] 0.725532165 0.462540989 0.033006719
> pf(x, 3, 1e6)
[1] 0.18784423 0.60837436 0.97585435
> pf(x, 3, Inf)
[1] 0.18784423 0.60837482 0.97585479
> 
> df(x, 1e6, 5)
[1] 0.158602071 0.610206081 0.061036395
> df(x, Inf, 5)
[1] 0.15859792 0.61020761 0.06103637
> pf(x, 1e6, 5)
[1] 0.0077295711 0.4158807972 0.9022692409
> pf(x, Inf, 5)
[1] 0.0077292503 0.4158801870 0.9022693759
> 
> df(x, Inf, Inf)# (0, Inf, 0)  - since 2.1.1
[1]   0 Inf   0
> pf(x, Inf, Inf)# (0, 1/2, 1)
[1] 0.0 0.5 1.0
> 
> pf(x, 5, Inf, ncp=0)
[1] 0.097730624 0.584119813 0.992270750
> all.equal(pf(x, 5, 1e6, ncp=1), tolerance = 1e-6,
+           c(0.065933194, 0.470879987, 0.978875867))
[1] TRUE
> all.equal(pf(x, 5, 1e7, ncp=1), tolerance = 1e-6,
+           c(0.06593309, 0.47088028, 0.97887641))
[1] TRUE
> all.equal(pf(x, 5, 1e8, ncp=1), tolerance = 1e-6,
+           c(0.0659330751, 0.4708802996, 0.9788764591))
[1] TRUE
> pf(x, 5, Inf, ncp=1)
[1] 0.065933078 0.470880318 0.978876467
> 
> dt(1, Inf)
[1] 0.24197072
> dt(1, Inf, ncp=0)
[1] 0.24197072
> dt(1, Inf, ncp=1)
[1] 0.39894228
> dt(1, 1e6, ncp=1)
[1] 0.39894208
> dt(1, 1e7, ncp=1)
[1] 0.39894226
> dt(1, 1e8, ncp=1)
[1] 0.39894227
> dt(1, 1e10, ncp=1) # = Inf
[1] 0.39894228
> ## Inf valid as from 2.1.1: df(x, 1e16, 5) was way off in 2.0.1.
> 
> sml.x <- c(10^-c(2:8,100), 0)
> cbind(x = sml.x, `dt(x,*)` = dt(sml.x, df = 2, ncp=1))
           x    dt(x,*)
 [1,]  1e-02 0.21686052
 [2,]  1e-03 0.21468294
 [3,]  1e-04 0.21446517
 [4,]  1e-05 0.21444339
 [5,]  1e-06 0.21444121
 [6,]  1e-07 0.21444100
 [7,]  1e-08 0.21444097
 [8,] 1e-100 0.21444097
 [9,]  0e+00 0.21444097
> ## small 'x' used to suffer from cancellation
> options(oo)
> x <- c(outer(1:12, 10^c(-3:2, 6:9, 10*(2:30))))
> for(nu in c(.75, 1.2, 4.5, 999, 1e50)) {
+     lfx <- dt(x, df=nu, log=TRUE)
+     stopifnot(is.finite(lfx), All.eq(exp(lfx), dt(x, df=nu)))
+ }## dt(1e160, 1.2, log=TRUE) was -Inf  up to R 2.15.2
> 
> ## pf() with large df1 or df2
> ## (was said to be PR#7099, but that is about non-central pchisq)
> nu <- 2^seq(25, 34, 0.5)
> target <- pchisq(1, 1) # 0.682...
> y <- pf(1, 1, nu)
> stopifnot(All.eq(pf(1, 1, Inf), target),
+           diff(c(y, target)) > 0, # i.e. pf(1, 1, *) is monotone increasing
+           abs(y[1] - (target - 7.21129e-9)) < 1e-11) # computed value
> ## non-monotone in R <= 2.1.0
> 
> stopifnot(pgamma(Inf, 1.1) == 1)
> ## didn't not terminate in R 2.1.x (only)
> 
> ## qgamma(q, *) should give {0,Inf} for q={0,1}
> sh <- c(1.1, 0.5, 0.2, 0.15, 1e-2, 1e-10)
> stopifnot(Inf == qgamma(1, sh))
> stopifnot(0   == qgamma(0, sh))
> ## the first gave Inf, NaN, and 99.425 in R 2.1.1 and earlier
> 
> ## In extreme left tail {PR#11030}
> p <- 10:123*1e-12
> qg <- qgamma(p, shape=19)
> qg2<- qgamma(1:100 * 1e-9, shape=11)
> stopifnot(diff(qg, diff=2) < -6e-6,
+           diff(qg2,diff=2) < -6e-6,
+ 	  abs(1 - pgamma(qg, 19)/ p) < 1e-13,
+           All.eq(qg  [1], 2.35047385139143),
+           All.eq(qg2[30], 1.11512318734547))
> ## was non-continuous in R 2.6.2 and earlier
> 
> f2 <- c(0.5, 1:4)
> stopifnot(df(0, 1, f2) == Inf,
+           df(0, 2, f2) == 1,
+           df(0, 3, f2) == 0)
> ## only the last one was ok in R 2.2.1 and earlier
> 
> x0 <- -2 * 10^-c(22,10,7,5) # ==> d*() warns about non-integer:
> assertWarning(fx0 <- dbinom(x0, size = 3, prob = 0.1))
> stopifnot(fx0 == 0,  pbinom(x0, size = 3, prob = 0.1) == 0)
> 
> ## very small negatives were rounded to 0 in R 2.2.1 and earlier
> 
> ## dbeta(*, ncp):
> db.x <- c(0, 5, 80, 405, 1280, 3125, 6480, 12005, 20480, 32805,
+ 	  50000, 73205, 103680, 142805, 192080, 253125, 327680)
> a <- rlnorm(100)
> stopifnot(All.eq(a, dbeta(0, 1, a, ncp=0)),
+ 	  dbeta(0, 0.9, 2.2, ncp = c(0, a)) == Inf,
+ 	  All.eq(65536 * dbeta(0:16/16, 5,1), db.x),
+ 	  All.eq(exp(16 * log(2) + dbeta(0:16/16, 5,1, log=TRUE)), db.x)
+           )
> ## the first gave 0, the 2nd NaN in R <= 2.3.0; others use 'TRUE' values
> stopifnot(all.equal(dbeta(0.8, 0.5, 5, ncp=1000),# was way too small in R <= 2.6.2
+ 		    3.001852308909e-35),
+ 	  all.equal(1, integrate(dbeta, 0,1, 0.8, 0.5, ncp=1000)$value,
+ 		    tolerance = 1e-4),
+           all.equal(1, integrate(dbeta, 0,1, 0.5, 200, ncp=720)$value),
+           all.equal(1, integrate(dbeta, 0,1, 125, 200, ncp=2000)$value)
+           )
> 
> ## df(*, ncp):
> x <- seq(0, 10, length=101)
> h <- 1e-7
> dx.h <- (pf(x+h, 7, 5, ncp= 2.5) - pf(x-h, 7, 5, ncp= 2.5)) / (2*h)
> stopifnot(all.equal(dx.h, df(x, 7, 5, ncp= 2.5), tolerance = 1e-6),# (1.50 | 1.65)e-8
+           All.eq(df(0, 2, 4, ncp=x), df(1e-300, 2, 4, ncp=x))
+           )
> 
> ## qt(p ~ 0, df=1) - PR#9804
> p <- 10^(-10:-20)
> qtp <- qt(p, df = 1)
> ## relative error < 10^-14 :
> stopifnot(abs(1 - p / pt(qtp, df=1)) < 1e-14)
> 
> ## Similarly for df = 2 --- both for p ~ 0  *and*  p ~ 1/2
> ## P ~ 0
> stopifnot(all.equal(qt(-740, df=2, log=TRUE), -exp(370)/sqrt(2)))
> ## P ~ 1 (=> p ~ 0.5):
> p.5 <- 0.5 + 2^(-5*(5:8))
> stopifnot(all.equal(qt(p.5, df = 2),
+ 		    c(8.429369702179e-08, 2.634178031931e-09,
+ 		      8.231806349784e-11, 2.572439484308e-12)))
> ## qt(<large>, log = TRUE)  is now more finite and monotone (again!):
> stopifnot(all.equal(qt(-1000, df = 4, log=TRUE),
+                     -4.930611e108, tolerance = 1e-6))
> qtp <- qt(-(20:850), df=1.2, log=TRUE, lower=FALSE)
> ##almost: stopifnot(all(abs(5/6 - diff(log(qtp))) < 1e-11))
> stopifnot(abs(5/6 - quantile(diff(log(qtp)), pr=c(0,0.995))) < 1e-11)
> 
> ## close to df=1 (where Taylor steps are important!):
> stopifnot(all.equal(-20, pt(qt(-20, df=1.02, log=TRUE),
+                           df=1.02, log=TRUE), tolerance = 1e-12),
+           diff(lq <- log(qt(-2^-(10:600), df=1.1, log=TRUE))) > 0.6)
> lq1 <- log(qt(-2^-(20:600), df=1, log=TRUE))
> lq2 <- log(qt(-2^-(20:600), df=2, log=TRUE))
> stopifnot(mean(abs(diff(lq1) - log(2)      )) < 1e-8,
+ 	  mean(abs(diff(lq2) - log(sqrt(2)))) < 4e-8)
> ## Case, where log.p=TRUE was fine, but log.p=FALSE (default) gave NaN:
> lp <- 40:406
> stopifnot(all.equal(lp, -pt(qt(exp(-lp), 1.2), 1.2, log=TRUE), tolerance = 4e-16))
> 
> 
> ## pbeta(*, log=TRUE) {toms708} -- now improved tail behavior
> x <- c(.01, .10, .25, .40, .55, .71, .98)
> pbval <- c(-0.04605755624088, -0.3182809860569, -0.7503593555585,
+            -1.241555830932, -1.851527837938, -2.76044482378, -8.149862739881)
> stopifnot(all.equal(pbeta(x, 0.8, 2, lower=FALSE, log=TRUE), pbval),
+           all.equal(pbeta(1-x, 2, 0.8, log=TRUE), pbval))
> qq <- 2^(0:1022)
> df.set <- c(0.1, 0.2, 0.5, 1, 1.2, 2.2, 5, 10, 20, 50, 100, 500)
> for(nu in df.set) {
+     pqq <- pt(-qq, df = nu, log=TRUE)
+     stopifnot(is.finite(pqq))
+ }
> ## PR#14230 -- more extreme beta cases {should no longer rely on denormalized}
> x <- (256:512)/1024
> P <- pbeta(x, 3, 2200, lower.tail=FALSE, log.p=TRUE)
> stopifnot(is.finite(P), P < -600,
+ 	  -.001 < (D3P <- diff(P, diff = 3)), D3P < 0, diff(D3P) < 0)
> ## all but the first 43 where -Inf in R <= 2.9.1
> stopifnot(All.eq(pt(2^-30, df=10),
+                  0.50000000036238542))
> ## = .5+ integrate(dt, 0,2^-30, df=10, rel.tol=1e-20)
> 
> ## rbinom(*, size) gave NaN for large size up to R <= 2.6.1
> M <- .Machine$integer.max
> set.seed(7)
> tt <- table(rbinom(100,    M, pr = 1e-9)) # had values in {0,2} only
> t2 <- table(rbinom(100, 10*M, pr = 1e-10))
> stopifnot(names(tt) == 0:6, sum(tt) == 100, sum(t2) == 100) ## no NaN there
> 
> ## qf() with large df1, df2  and/or  small p:
> x <- 0.01; f1 <- 1e60; f2 <- 1e90
> stopifnot(qf(1/4, Inf, Inf) == 1,
+ 	  all.equal(1, 1e-18/ pf(qf(1e-18, 12,50), 12,50), tolerance = 1e-10),
+ 	  abs(x - qf(pf(x, f1,f2, log.p=TRUE), f1,f2, log.p=TRUE)) < 1e-4)
> 
> ## qbeta(*, log.p) for "border" case:
> stopifnot(is.finite(q0 <- qbeta(-1e10, 50,40, log.p=TRUE)),
+           1 ==            qbeta(-1e10,  2, 3, log.p=TRUE, lower=FALSE))
> ## infinite loop or NaN in R <= 2.7.0
> 
> ## phyper(x, 0,0,0), notably for huge x
> stopifnot(all(phyper(c(0:3, 1e67), 0,0,0) == 1))
> ## practically infinite loop and NaN in R <= 2.7.1 (PR#11813)
> 
> ## plnorm(<= 0, . , log.p=TRUE)
> stopifnot(plnorm(-1:0, lower.tail=FALSE, log.p=TRUE) == 0,
+ 	  plnorm(-1:0, lower.tail=TRUE,	 log.p=TRUE) == -Inf)
> ## was wrongly == 'log.p=FALSE' up to R <= 2.7.1 (PR#11867)
> 
> 
> ## pchisq(df=0) was wrong in 2.7.1; then, upto 2.10.1, P*(0,0) gave 1
> stopifnot(pchisq(c(-1,0,1), df=0) == c(0,0,1),
+           pchisq(c(-1,0,1), df=0, lower.tail=FALSE) == c(1,1,0),
+ 	  ## for ncp >= 80, gave values >= 1 in 2.10.0
+ 	  pchisq(500:700, 1.01, ncp = 80) <= 1)
> 
> ## dnbinom for extreme  size and/or mu :
> mu <- 20
> d <- dnbinom(17, mu=mu, size = 1e11*2^(1:10)) - dpois(17, lambda=mu)
> stopifnot(d < 0, diff(d) > 0, d[1] < 1e-10)
> ## was wrong up to 2.7.1
> ## The fix to the above, for x = 0, had a new cancellation problem
> mu <- 1e12 * 2^(0:20)
> stopifnot(all.equal(1/(1+mu), dnbinom(0, size = 1, mu = mu), tolerance = 1e-13))
> ## was wrong in 2.7.2 (only)
> mu <- sort(outer(1:7, 10^c(0:10,50*(1:6))))
> NB <- dnbinom(5, size=1e305, mu=mu, log=TRUE)
> P  <- dpois  (5,                mu, log=TRUE)
> stopifnot(abs(rErr(NB,P)) < 9*Meps)# seen 2.5*
> ## wrong in 3.1.0 and earlier
> 
> 
> ## Non-central F for large x
> x <- 1e16 * 1.1 ^ (0:20)
> dP <- diff(pf(x, df1=1, df2=1, ncp=20, lower.tail=FALSE, log=TRUE))
> stopifnot(-0.047 < dP, dP < -0.0455)
> ## pf(*, log) jumped to -Inf prematurely in 2.8.0 and earlier
> 
> 
> ## Non-central Chi^2 density for large x
> stopifnot(0 == dchisq(c(Inf, 1e80, 1e50, 1e40), df=10, ncp=1))
> ## did hang in 2.8.0 and earlier (PR#13309).
> 
> 
> ## qbinom() .. particularly for large sizes, small prob:
> p.s <- c(.01, .001, .1, .25)
> pr <- (2:20)*1e-7
> sizes <- 1000*(5000 + c(0,6,16)) + 279
> k.s <- 0:15; q.xct <- rep(k.s, each=length(pr))
> for(sz in sizes) {
+     for(p in p.s) {
+ 	qb <- qbinom(p=p, size = sz, prob=pr)
+ 	pb <- qpois (p=p, lambda = sz * pr)
+ 	stopifnot(All.eq(qb, pb))
+     }
+     pp.x <- outer(pr, k.s, function(pr, q) pbinom(q, size = sz, prob=pr))
+     qq.x <- apply(pp.x, 2, function(p)     qbinom(p, size = sz, prob=pr))
+     stopifnot(qq.x == q.xct)
+ }
> ## do_search() in qbinom() contained a thinko up to 2.9.0 (PR#13711)
> 
> 
> ## pbeta(x, a,b, log=TRUE)  for small x and a  is ~ log-linear
> x <- 2^-(200:10)
> for(a in c(1e-8, 1e-12, 16e-16, 4e-16))
+     for(b in c(0.6, 1, 2, 10)) {
+         dp <- diff(pbeta(x, a, b, log=TRUE)) # constant approximately
+         stopifnot(sd(dp) / mean(dp) < 0.0007)
+     }
> ## had  accidental cancellation '1 - w'
> 
> ## qgamma(p, a) for small a and (hence) small p
> ## pgamma(x, a) for very very small a
> a <- 2^-seq(10,1000, .25)
> q.1c <- qgamma(1e-100,a,lower.tail=FALSE)
> q.3c <- qgamma(1e-300,a,lower.tail=FALSE)
> p.1c <- pgamma(q.1c[q.1c > 0], a[q.1c > 0], lower.tail=FALSE)
> p.3c <- pgamma(q.3c[q.3c > 0], a[q.3c > 0], lower.tail=FALSE)
> x <- 1+1e-7*c(-1,1); pg <- pgamma(x, shape = 2^-64, lower.tail=FALSE)
> stopifnot(qgamma(.99, .00001) == 0,
+           abs(pg[2] - 1.18928249197237758088243e-20) < 1e-33,
+ 	  abs(diff(pg) + diff(x)*dgamma(1, 2^-64)) < 1e-13 * mean(pg),
+ 	  abs(1 - p.1c/1e-100) < 10e-13,# max = 2.243e-13 / 2.442 e-13
+ 	  abs(1 - p.3c/1e-300) < 28e-13)# max = 7.057e-13
> ## qgamma() was wrong here, orders of magnitude up to R 2.10.0
> ## pgamma() had inaccuracies, e.g.,
> ## pgamma(x, shape = 2^-64, lower.tail=FALSE)  was discontinuous at x=1
> 
> stopifnot(all(qpois((0:8)/8, lambda=0) == 0))
> ## gave Inf as p==1 was checked *before* lambda==0
> 
> ## extreme tail of non-central chisquare
> stopifnot(all.equal(pchisq(200, 4, ncp=.001, log.p=TRUE), -3.851e-42))
> ## jumped to zero too early up to R 2.10.1 (PR#14216)
> ## left "extreme tail"
> lp <- pchisq(2^-(0:200), 100, 1, log=TRUE)
> stopifnot(is.finite(lp), lp < -184,
+ 	  all.equal(lp[201], -7115.10693158))
> dlp <- diff(lp)
> dd <- abs(dlp[-(1:30)] - -34.65735902799)
> stopifnot(-34.66 < dlp, dlp < -34.41, dd < 1e-8)# 2.2e-10 64bit Lnx
> ## underflowed to -Inf much too early in R <= 3.1.0
> for(e in c(0, 2e-16))# continuity at 80 (= branch point)
+ stopifnot(all.equal(pchisq(1:2, 1.01, ncp = 80*(1-e), log=TRUE),
+ 		    c(-34.57369629, -31.31514671)))
> 
> ## logit() == qlogit() on the right extreme:
> x <- c(10:80, 80 + 5*(1:24), 200 + 20*(1:25))
> stopifnot(All.eq(x, qlogis(plogis(x, log.p=TRUE),
+ 			   log.p=TRUE)))
> ## qlogis() gave Inf much too early for R <= 2.12.1
> ## Part 2:
> x <- c(x, seq(700, 800, by=10))
> stopifnot(All.eq(x, qlogis(plogis(x, lower=FALSE, log.p=TRUE),
+ 			   lower=FALSE, log.p=TRUE)))
> # plogis() underflowed to -Inf too early for R <= 2.15.0
> 
> ## log upper tail pbeta():
> x <- (25:50)/128
> pbx <- pbeta(x, 1/2, 2200, lower.tail=FALSE, log.p=TRUE)
> d2p <- diff(dp <- diff(pbx))
> b <- 2200*2^(0:50)
> y <- log(-pbeta(.28, 1/2, b, lower.tail=FALSE, log.p=TRUE))
> stopifnot(-1094 < pbx, pbx < -481.66,
+             -29 < dp,   dp < -20,
+            -.36 < d2p, d2p < -.2,
+           all.equal(log(b), y+1.113, tolerance = .00002)
+           )
> ## pbx had two -Inf; y was all Inf  for R <= 2.15.3;  PR#15162
> 
> ## dnorm(x) for "large" |x|
> stopifnot(abs(1 - dnorm(35+3^-9)/ 3.933395747534971e-267) < 1e-15)
> ## has been losing up to 8 bit precision for R <= 3.0.x
> 
> ## pbeta(x, <small a>,<small b>, .., log):
> ldp <- diff(log(diff(pbeta(0.5, 2^-(90+ 1:25), 2^-60, log.p=TRUE))))
> stopifnot(abs(ldp - log(1/2)) < 1e-9)
> ## pbeta(*, log) lost all precision here, for R <= 3.0.x (PR#15641)
> ##
> ## "stair function" effect (from denormalized numbers)
> a <- 43779; b <- 0.06728
> x. <- .9833 + (0:100)*1e-6
> px <- pbeta(x., a,b, log=TRUE) # plot(x., px) # -> "stair"
> d2. <- diff(dpx <- diff(px))
> stopifnot(all.equal(px[1], -746.0986886924, tol=1e-12),
+           0.0445741 < dpx, dpx < 0.0445783,
+           -4.2e-8 < d2., d2. < -4.18e-8)
> ## were way off in R <= 3.1.0
> 
> c0 <- system.time(p0 <- pbeta(  .9999, 1e30, 1.001, log=TRUE))
> cB <- max(.001, c0[[1]])# base time
> c1 <- system.time(p1 <- pbeta(1- 1e-9, 1e30, 1.001, log=TRUE))
> c2 <- system.time(p2 <- pbeta(1-1e-12, 1e30, 1.001, log=TRUE))
> stopifnot(all.equal(p0, -1.000050003333e26, tol=1e-10),
+           all.equal(p1, -1e21, tol = 1e-6),
+           all.equal(p2, -9.9997788e17),
+           c(c1[[1]], c2[[1]]) < 1000*cB)
> ## (almost?) infinite loop in R <= 3.1.0
> 
> 
> ## pbinom(), dbinom(), dhyper(),.. : R allows "almost integer" n
> for (FUN in c(function(n) dbinom(1,n,0.5), function(n) pbinom(1,n,0.5),
+               function(n) dpois(n, n), function(n) dhyper(n+1, n+5,n+5, n)))
+     try( lapply(sample(10000, size=1000), function(M) {
+     ## invisible(lapply(sample(10000, size=1000), function(M) {
+         n <- (M/100)*10^(2:20); if(anyNA(P <- FUN(n)))
+             stop("NA for M=",M, "; 10ex=",paste((2:20)[is.na(P)], collapse=", "))}))
> ## check was too tight for large n in R <= 3.1.0 (PR#15734)
> 
> ## [dpqr]beta(*, a,b) where a and/or b are Inf
> stopifnot(pbeta(.1, Inf, 40) == 0,
+           pbeta(.5, 40, Inf) == 1,
+           pbeta(.4, Inf,Inf) == 0,
+           pbeta(.5, Inf,Inf) == 1,
+           ## gave infinite loop (or NaN) in R <= 3.1.0
+           qbeta(.9, Inf, 100) == 1, # Inf.loop
+           qbeta(.1, Inf, Inf) == 1/2)# NaN + Warning
> ## range check (in "close" cases):
> assertWarning(qN <- qbeta(2^-(10^(1:3)), 2,3, log.p=TRUE))
> assertWarning(qn <- qbeta(c(-.1, -1e-300, 1.25), 2,3))
> stopifnot(is.nan(qN), is.nan(qn))
> 
> ## lognormal boundary case sdlog = 0:
> p <- (0:8)/8; x <- 2^(-10:10)
> stopifnot(all.equal(qlnorm(p, meanlog=1:2, sdlog=0),
+ 		    qlnorm(p, meanlog=1:2, sdlog=1e-200)),
+ 	  dlnorm(x, sdlog=0) == ifelse(x == 1, Inf, 0))
> 
> ## qbeta(*, a,b) when  a,b << 1 : can easily fail
> qbeta(2^-28, 0.125, 2^-26) # 1000 Newton it + warning
[1] 9.094955e-13
> a <- 1/8; b <- 2^-(4:200); alpha <- b/4
> qq <- qbeta(alpha, a,b)# gave warnings intermediately
> pp <- pbeta(qq, a,b)
> stopifnot(pp > 0, diff(pp) < 0, ## pbeta(qbeta(alpha,*),*) == alpha:
+           abs(1 - pp/alpha) < 4e-15)# seeing 2.2e-16
> 
> ## orig. qbeta() using *many* Newton steps; case where we "know the truth"
> a <- 25; b <- 6; x <- 2^-c(3:15, 100, 200, 250, 300+100*(0:7))
> pb <- c(## via Rmpfr's roundMpfr(pbetaI(x, a,b, log.p=TRUE, precBits = 2048), 64) :
+     -40.7588797271766572448, -57.7574063441183625303, -74.9287878018119846216,
+     -92.1806244636893542185, -109.471318248524419364, -126.781111923947395655,
+     -144.100375042814531426, -161.424352961544612370, -178.750683324909148353,
+     -196.078188674895169383, -213.406281209657976525, -230.734667259724367416,
+     -248.063200048177428608, -1721.00081201679567511, -3453.86876341665894863,
+     -4320.30273911659058550, -5186.73671481652222237, -6919.60466621638549567,
+     -8652.47261761624876897, -10385.3405690161120427, -12118.2085204159753165,
+     -13851.0764718158385902, -15583.9444232157018631, -17316.8123746155651368)
> stopifnot(all.equal(pb, pbeta(x,a,b, log.p=TRUE), tol=8e-16))# seeing {1.5|1.6|2.0}e-16
> qp <- qbeta(pb, a,b, log.p=TRUE)
> ## x == qbeta(pbeta(x, *), *) :
> stopifnot(qp > 0, all.equal(x, qp, tol= 1e-15))# seeing {2.4|3.3}e-16
> 
> ## qbeta(), PR#15755
> a1 <- 0.0672788; b1 <- 226390
> p <- 0.6948886
> qp <- qbeta(p, a1,b1)
> stopifnot(qp < 2e-8, # was '1' (with a warning) in R <= 3.1.0
+           All.eq(p, pbeta(qp, a1,b1)))
> ## less extreme example, same phenomenon:
> a <- 43779; b <- 0.06728
> stopifnot(All.eq(0.695, pbeta(qbeta(0.695, b,a), b,a)))
> x <- -exp(seq(0, 14, by=2^-9))
> ct <- system.time(qx <- qbeta(x, a,b, log.p=TRUE))[[1]]
> pqx <- pbeta(qx, a,b, log=TRUE)
> stopifnot(all.equal(x, pqx, tol= 2e-15)) # seeing {3.51|3.54}e-16
> ## note that qx[x > -exp(2)] is too close to 1 to get full accuracy:
> ## i2 <- x > -exp(2); all.equal(x[i2], pqx[i2], tol= 0)#-> 5.849e-12
> if(ct > 0.5) { cat("system.time:\n"); print(ct) }# lynne(2013): 0.048
> ## was Inf, and much slower, for R <= 3.1.0
> x3 <- -(15450:15700)/2^11
> pq3 <- pbeta(qbeta(x3, a,b, log.p=TRUE), a,b, log=TRUE)
> stopifnot(mean(abs(pq3-x3)) < 4e-12,# 1.46e-12
+           max (abs(pq3-x3)) < 8e-12)# 2.95e-12
> ##
> .a <- .2; .b <- .03; lp <- -(10^-(1:323))
> qq <- qbeta(lp, .a,.b, log=TRUE) # warnings in R <= 3.1.0
> assertWarning(qN <- qbeta(.5, 2,3, log.p=TRUE))
> assertWarning(qn <- qbeta(c(-.1, 1.25), 2,3))
> stopifnot(1-qq < 1e-15, is.nan(qN), is.nan(qn))# typically qq == 1  exactly
> ## failed in intermediate versions
> ##
> a <- 2^-8; b <- 2^(200:500)
> pq <- pbeta(qbeta(1/8, a, b), a, b)
> stopifnot(abs(pq - 1/8) < 1/8)
> ## whereas  qbeta() would underflow to 0 "too early", for R <= 3.1.0
> #
> ## very extreme tails on both sides
> x <- c(1e-300, 1e-12, 1e-5, 0.1, 0.21, 0.3)
> stopifnot(0 == qbeta(x, 2^-12, 2^-10))## gave warnings
> a <- 10^-(8:323)
> qb <- qbeta(0.95, a, 20)
> ## had warnings and wrong value +1; also NaN
> ct2 <- system.time(q2 <- qbeta(0.95, a,a))[1]
> stopifnot(is.finite(qb), qb < 1e-300, q2 == 1)
> if(ct2 > 0.020) { cat("system.time:\n"); print(ct2) }
> ## had warnings and was much slower for R <= 3.1.0
> 
> ## qt(p, df= Inf, ncp)  <==> qnorm(p, m=ncp)
> p <- (0:32)/32
> stopifnot(all.equal(qt(p, df=Inf, ncp=5), qnorm(p, m=5)))
> ## qt(*, df=Inf, .)  gave NaN in  R <= 3.2.1
> 
> ## rhyper(*, <large>);  PR#16489
> ct3 <- system.time(N <- rhyper(100, 8000, 1e9-8000, 1e6))[1]
> table(N)
N
 0  2  3  4  5  6  7  8  9 10 11 12 13 14 15 
 1  1  2  2  8 21  9 18 15 10  1  7  1  2  2 
> summary(N)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    0.0     6.0     8.0     7.9     9.0    15.0 
> stopifnot(abs(mean(N) - 8) < 1.5)
> if(ct3 > 0.02) { cat("system.time:\n"); print(ct3) }
> ## N were all 0 and took very long for R <= 3.2.1
> set.seed(17)
> stopifnot(rhyper(1, 3024, 27466, 251) == 25,
+           rhyper(1,  329,  3059, 225) == 22)
> ## failed for a day after a "thinko" in the above bug fix.
> 
> ## *chisq(*, df=0, ncp=0) == Point mass at 0
> stopifnot(rchisq(32,        df=0, ncp=0) == 0,
+           dchisq((0:16)/16, df=0, ncp=0) == c(Inf, rep(0, 16)))
> ## gave all NaN's  for R <= 3.2.2
> 
> ## pchisq(*, df=0, ncp > 0, log.p=TRUE) :
> th <- 10*c(1:9,10^c(1:3,7))
> pp <- pchisq(0, df = 0, ncp=th, log.p=TRUE)
> stopifnot(all.equal(pp, -th/2, tol=1e-15))
> ## underflowed at about th ~= 60  in R <= 3.2.2
> 
> ## pnbinom (-> C's bratio())
> op <- options(warn = 1)# -- NaN's giving warnings
> L <- 1e308; p <- suppressWarnings(pnbinom(L, L, mu = 5)) # NaN or 1 (for 64 / 32 bit)
> is.nan(p) || p == 1
[1] TRUE
> ## gave infinite loop on some 64b platforms in R <= 3.2.3
> 
> ## [dpqr]nbinom(*, mu, size=Inf) -- PR#16727
> L <- 1e308; mu <- 5; pp <- (0:16)/16
> x <- c(0:3, 1e10, 1e100, L, Inf)
> (d <- dnbinom(x,  mu = mu, size = Inf)) # gave NaN (for 0 and L)
[1] 0.006737947 0.033689735 0.084224337 0.140373896 0.000000000 0.000000000
[7] 0.000000000 0.000000000
> (p <- pnbinom(x,  mu = mu, size = Inf)) # gave all NaN
[1] 0.006737947 0.040427682 0.124652019 0.265025915 1.000000000 1.000000000
[7] 1.000000000 1.000000000
> (q <- qnbinom(pp, mu = mu, size = Inf)) # gave all NaN
 [1]   0   2   3   3   3   4   4   4   5   5   6   6   6   7   8   9 Inf
> set.seed(1); NI <- rnbinom(32, mu = mu, size = Inf)# gave all NaN
> set.seed(1); N2 <- rnbinom(32, mu = mu, size = L  )
> stopifnot(all.equal(d, dpois(x, mu)),
+ 	  all.equal(p, ppois(x, mu)),
+ 	  q == qpois(pp, mu),
+ 	  identical(NI, N2))
> options(op)
> ## size = Inf -- mostly gave NaN  in R <= 3.2.3
> 
> ## qpois(p, *) for invalid 'p' should give NaN -- PR#16972
> stopifnot(is.nan(suppressWarnings(c(qpois(c(-2,3, NaN), 3), qpois(1, 3, log.p=TRUE),
+                                     qpois(.5, 0, log.p=TRUE), qpois(c(-1,pi), 0)))))
> ## those in the 2nd line gave 0 in R <= 3.3.1
> ## Similar but different for qgeom():
> stopifnot(qgeom((0:8)/8, prob=1) == 0, ## p=1 gave Inf in R <= 3.3.1
+           is.nan(suppressWarnings(qgeom(c(-1/4, 1.1), prob=1))))
> 
> ## all our RNG  r<dist>() functions:
> ##' catch all: value and warnings or error <-- demo(error.catching) :
> tryCatch.W.E <- function(expr) {
+     W <- NULL
+     w.handler <- function(w){ # warning handler
+ 	W <<- w
+ 	invokeRestart("muffleWarning")
+     }
+     list(value = withCallingHandlers(tryCatch(expr, error = function(e) e),
+ 				     warning = w.handler),
+ 	 warning = W)
+ }
> .stat.ns <- asNamespace("stats")
> Ns <- 4
> for(dist in PDQR) {
+     fn <- paste0("r",dist)
+     cat(sprintf("%-9s(%d, ..): ", fn, Ns))
+     F <- get(fn, envir = .stat.ns)
+     nArg <- length(fms <- formals(F))
+     if(dist %in% c("nbinom", "gamma")) ## cannot specify *both* 'prob' & 'mu' / 'rate' & 'scale'
+         nArg <- nArg - 1
+     nA1 <- nArg - 1 # those beside the first (= 'n' mostly)
+     expected <- rep(if(dist %in% PDQRinteg) NA_integer_ else NaN, Ns)
+     for(ia in seq_len(nA1)) {
+         aa <- rep(list(1), nA1)
+         aa[[ia]] <- NA
+         cat(ia,"")
+         R <- tryCatch.W.E( do.call(F, c(Ns, aa)) )
+         if(!inherits(R$warning, "simpleWarning")) cat(" .. did *NOT* give a warning! ")
+ 	if(!(identical(R$value, expected))) { ## allow NA/NaN mismatch in these cases for now:
+ 	    if(!(dist %in% c("beta","f","t") && all(is.na(R$value))))
+ 		cat(" .. not giving expected NA/NaN's ")
+         }
+     }
+     cat(" [Ok]\n")
+ }
rbinom   (4, ..): 1 2  [Ok]
rgeom    (4, ..): 1  [Ok]
rhyper   (4, ..): 1 2 3  [Ok]
rnbinom  (4, ..): 1 2  [Ok]
rpois    (4, ..): 1  [Ok]
rsignrank(4, ..): 1  [Ok]
rwilcox  (4, ..): 1 2  [Ok]
rbeta    (4, ..): 1 2 3  [Ok]
rcauchy  (4, ..): 1 2  [Ok]
rchisq   (4, ..): 1 2  [Ok]
rexp     (4, ..): 1  [Ok]
rf       (4, ..): 1 2 3  [Ok]
rgamma   (4, ..): 1 2  [Ok]
rlnorm   (4, ..): 1 2  [Ok]
rlogis   (4, ..): 1 2  [Ok]
rnorm    (4, ..): 1 2  [Ok]
rt       (4, ..): 1 2  [Ok]
runif    (4, ..): 1 2  [Ok]
rweibull (4, ..): 1 2  [Ok]
> 
> 
> ## qbeta() in very asymmetric cases
> sh2 <- 2^seq(9,16, by=1/16)
> qbet <- qbeta(1e-10, 1.5, shape2=sh2, lower.tail=FALSE)
> plot(sh2, 1- pbeta(qbet, 1.5, sh2, lower.tail=FALSE) * 1e10, log="x")
> dqb <- diff(qbet); d2qb <- diff(dqb); d3qb <- diff(d2qb)
> stopifnot(all.equal(qbet[[1]], 0.047206901483498, tol=1e-12),
+ 	  max(abs(1- pbeta(qbet, 1.5, sh2, lower.tail=FALSE) * 1e10)) < 1e-12,# Lx 64b: 2.4e-13
+ 	  0 > dqb, dqb > -0.002,
+ 	  0 < d2qb, d2qb < 0.00427,
+ 	  -3.2e-8 > d3qb, d3qb > -3.1e-6,
+ 	  diff(d3qb) > 1e-9)
> ## had discontinuity (from wrong jump out of Newton) in R <= 3.3.2
> 
> 
> ## rt() [PR#17306];  rf() and rbeta() [PR#17375] with non-scalar 'ncp'
> nc <- c(NA, 1); iN <- is.na(rep_len(nc, 3))
> ## each gives warning  "NAs produced":
> assertWarning(T <- rt   (3, 4,   ncp = nc))
> assertWarning(F <- rf   (3, 4,5, ncp = nc))
> assertWarning(B <- rbeta(3, 4,5, ncp = nc))
> stopifnot(identical(iN, is.na(T)), identical(iN, is.na(F)), identical(iN, is.na(B)))
> ## was not handled correctly, notably with NA's in ncp, in R <= 3.4.(2|3)
> 
> 
> ## check old version of walker_Probsample is being used for old sample kind
> suppressWarnings(RNGversion("3.5.0"))
> set.seed(12345)
> p <- c(2, rep(1, 200))
> x <- sample(length(p), 100000, prob = p, replace = TRUE)
> stopifnot(sum(x == 1) == 994)
> 
> ## check for faiure of new walker_Probsample
> RNGversion("3.6.0")
> set.seed(12345)
> epsilon <- 1e-10
> p201 <- prop.table( rep( c(1, epsilon), c(201, 999-201)))
> x <- sample(length(p201), 100000, prob = p201, replace = TRUE)
> stopifnot(sum(x <= 201) == 100000)
> 
> 
> cat("Time elapsed: ", proc.time() - .ptime,"\n")
Time elapsed:  1.522 0.014 1.542 0 0 
>