R version 2.4.0 Under development (unstable) (2006-05-20 r38139) Copyright (C) 2006 The R Foundation for Statistical Computing ISBN 3-900051-07-0 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 > > ###-- 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 = .Options$rErr.eps) + { + if(is.null(eps)) { eps <- 1e-30; options(rErr.eps = eps) } + 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()) + .Random.seed <- c(0,rep(7654, 3)) > > ## 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)) { + cat("n=",n,": ") + for(p in c(0,1,rbeta(n2, 2,4))) { + cat(".") + 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: + tst1 <- all.equal( pbinom(0:k, size = n, prob = p), + cumsum(dbinom(0:k, size = n, prob = p))) + tst <- 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))) + if(!isTRUE(tst1) || !isTRUE(tst)) { + cat("n=", n,"; p =",format(p),". k =",k) + if(!isTRUE(tst1)) cat("; tst1=",tst1) + if(!isTRUE(tst )) cat("; tst=", tst) + cat("\n") + } + } + } + cat("\n") + } n= 45 : ...................... n= 35 : ...................... n= 32 : ...................... n= 46 : ...................... n= 42 : ...................... n= 42 : ...................... n= 35 : ...................... n= 46 : ...................... n= 42 : ...................... n= 33 : ...................... n= 33 : ...................... n= 47 : ...................... n= 37 : ...................... n= 34 : ...................... n= 51 : ...................... n= 48 : ...................... > > ##__ 2. Geometric __ > for(pr in seq(1e-10,1,len=15)) { # p=0 is not a distribution + print(All.eq((dg <- dgeom(0:10, pr)), + pr * (1-pr)^(0:10))) + print(All.eq(cumsum(dg), pgeom(0:10, pr))) + } [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE > > ##__ 3. Hypergeometric __ > > m <- 10; n <- 7 > for(k in 2:m) { + x <- 0:(k+1) + print(All.eq(phyper(x, m, n, k), cumsum(dhyper(x, m, n, k)))) + } [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE > > ##__ 4. Negative Binomial __ > > ## PR #842 > for(size in seq(0.8,2, by=.1)) + print(all.equal(cumsum(dnbinom(0:7, size, .5)), + pnbinom(0:7, size, .5))) [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE > All.eq(pnbinom(c(1,3), .9, .5), c(0.777035760338812, 0.946945347071519)) [1] TRUE > > ##__ 5. Poisson __ > > all(dpois(0:5,0) == c(1, rep(0,5))) [1] TRUE > all(dpois(0:5,0, log=TRUE) == c(0, rep(-Inf, 5))) [1] TRUE > > ## 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)) { + tst <- all.equal(1 - pchisq(2*lambda, 2*(1+ 0:k)), + pp <- cumsum(dpois(0:k, lambda=lambda)), tol= 100*Meps) + if(!isTRUE(tst)) + cat("lambda=", format(lambda),". k =",k, " --> tst=", tst,"\n") + tst2 <- all.equal(pp, ppois(0:k, lambda=lambda), tol = 100*Meps) + if(!isTRUE(tst2)) + cat("lambda=", format(lambda),". k =",k, " --> tst2=", tst2,"\n") + tst3 <- all.equal(1 - pp, ppois(0:k, lambda=lambda, lower.tail=FALSE)) + if(!isTRUE(tst3)) + cat("lambda=", format(lambda),". k =",k, " --> tst3=", tst3,"\n") + } > > > ##__ 6. SignRank __ > for(n in rpois(32, lam=8)) { + x <- -1:(n + 4) + if(!isTRUE(eq <- All.eq(psignrank(x, n), cumsum(dsignrank(x, n))))) + print(eq) + } > > ##__ 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)) + if(!isTRUE(eq <- All.eq(Fx, cumsum(fx)))) + print(eq) + } > is.sym [1] TRUE > > > ###-------- 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)) { + tst <- all.equal((d1 <- dgamma( x, shape = sh, scale = sig)), + (d2 <- dgamma(x/sig, shape = sh, scale = 1) / sig), + tol = 1e-14)## __ad interim__ was 1e-15 + if(!isTRUE(tst)) + cat("ERROR: dgamma() doesn't scale:",tst,"\n", + " x =", formatC(x),"\n shape,scale=",formatC(c(sh, sig)),"\n") + tst <- All.eq(d1, (d3 <- 1/(Ga * sig^sh) * x^(sh-1) * exp(-x/sig))) + if(!isTRUE(tst)) + cat("NOT Equal:",tst,"\n x =", formatC(x), + "\n shape,scale=",formatC(c(sh, sig)),"\n") + } + } > pgamma(1,Inf,scale=Inf) == 0 [1] TRUE > ## Also pgamma(Inf,Inf) == 1 for which NaN was slightly more appropriate > all(is.nan(c(pgamma(Inf, 1,scale=Inf), + pgamma(Inf,Inf,scale=Inf)))) [1] TRUE Warning messages: 1: NaNs produced in: pgamma(q, shape, scale, lower.tail, log.p) 2: NaNs produced in: pgamma(q, shape, scale, lower.tail, log.p) > 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), tol=2e-15) + ) > > p <- 7e-4; df <- 0.9 > 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(log(1-p),df,log=T,lower=F),df, log=T,lower=F)/log(1-p) + )# 32b-i386: (2.2e-16, 0,0, 3.3e-16); Opteron: (2.2e-16, 0,0, 2.2e-15) + ) < 1e-14 [1] TRUE TRUE TRUE TRUE > > ##-- 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) > all.equal(qchisq(0.025,31,ncp=1,lower.tail=FALSE),# inf.loop PR#875 + 49.7766246561514, tol= 1e-11) [1] TRUE > for(df in c(0.1, 0.5, 1.5, 4.7, 10, 20,50,100)) { + cat("df =", formatC(df, wid=3)) + 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) + print(all.equal(xx, qchisq(pp, df=df, ncp=1), tol = dtol))# TRUE + ##or print(mapply(rErr, xx, qchisq(pp, df=df,ncp=1)), digits = 3) + } df = 0.1[1] TRUE df = 0.5[1] TRUE df = 1.5[1] TRUE df = 4.7[1] TRUE df = 10[1] TRUE df = 20[1] TRUE df = 50[1] TRUE df = 100[1] TRUE > > ## p ~= 1 (<==> 1-p ~= 0) -- gave infinite loop in R <= 1.8.1 -- PR#6421 > options(warn=-1) # ignore warnings from R's version of log1p > 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 > all.equal(q0[iO], q1[iO], 1e-5) [1] TRUE > all.equal(p0[iO], psml[iO]) [1] TRUE > options(warn=0) > > ##--- Beta (need more): > > ## big a & b (PR #643) > summary(a <- rlnorm(20, 5.5)) Min. 1st Qu. Median Mean 3rd Qu. Max. 27.65 105.20 190.50 310.50 322.00 1611.00 > summary(b <- rlnorm(20, 6.5)) Min. 1st Qu. Median Mean 3rd Qu. Max. 187.9 346.8 630.2 1281.0 978.0 9382.0 > pab <- expand.grid(seq(0,1,by=.1), a, b) > p <- pab[,1]; a <- pab[,2]; b <- pab[,3] > all.equal(dbeta(p,a,b), exp(pab <- dbeta(p,a,b, log = TRUE)), tol = 1e-11) [1] TRUE > sample(pab, 50) [1] -Inf -4.338148e+01 -Inf -3.321764e+02 -5.222040e+01 [6] -2.302052e+03 -6.743669e+03 -Inf -5.069544e+02 -Inf [11] -5.995978e+02 -7.505186e+00 -2.086599e+01 -Inf -3.287371e+02 [16] -4.161389e+02 -2.995488e+02 -6.933932e+01 -1.583166e+02 -1.557598e+01 [21] 9.825994e-01 -1.883283e+04 -6.418225e+02 1.741584e+00 -6.910231e+02 [26] -1.224343e+03 -1.974485e+01 -1.324279e+03 -2.013099e+03 -4.333847e+01 [31] -Inf -1.087940e+04 -9.095702e+00 -5.083759e+03 -Inf [36] -2.054621e+03 -1.309311e+02 -9.883952e+01 -1.685633e+02 -5.312428e+01 [41] -3.565202e+00 -1.790596e+02 -Inf -4.265242e+03 -Inf [46] -2.259480e+01 -9.688007e+02 -6.797130e+02 -5.032581e+02 -2.174204e+03 > > > ##--- Normal (& Lognormal) : > > qnorm(0) == -Inf && qnorm(-Inf, log = TRUE) == -Inf [1] TRUE > qnorm(1) == Inf && qnorm(0, log = TRUE) == Inf [1] TRUE > > is.nan(qnorm(1.1)) && + is.nan(qnorm(-.1)) # + warn [1] TRUE Warning messages: 1: NaNs produced in: qnorm(p, mean, sd, lower.tail, log.p) 2: NaNs produced in: qnorm(p, mean, sd, lower.tail, log.p) > > x <- c(-Inf, -1e100, 1:6, 1e200, Inf) > rbind(d.s0 =dnorm(x,3,s=0), p.s0 = pnorm(x,3,s=0), + d.sI =dnorm(x,3,s=Inf), p.sI = pnorm(x,3,s=Inf)) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] d.s0 0 0.0 0.0 0.0 Inf 0.0 0.0 0.0 0.0 0 p.s0 0 0.0 0.0 0.0 1.0 1.0 1.0 1.0 1.0 1 d.sI 0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0 p.sI 0 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 1 > > ## 3 Test data from Wichura (1988) : > all.equal(qnorm(c( 0.25, .001, 1e-20)), + c(-0.6744897501960817, -3.090232306167814, -9.262340089798408), + tol = 1e-15) [1] TRUE > # extreme tail -- available on log scale only: > all.equal(qnorm(-1e5, log = TRUE), -447.1974945) [1] TRUE > > z <- rnorm(1000); all.equal(pnorm(z), 1 - pnorm(-z), tol= 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) if(!isTRUE(all.equal(pt(z, df), 1 - pt(-z,df), tol= 1e-15))) + cat("ERROR -- df = ", df, "\n") > All.eq(pz <- pnorm(z), 1 - pnorm(z, lower=FALSE)) [1] TRUE > All.eq(pz, pnorm(-z, lower=FALSE)) [1] TRUE > All.eq(log(pz[z.ok]), pnorm(z[z.ok], log=TRUE)) [1] 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 > All.eq(pz, plnorm(exp(z))) [1] TRUE > > > ###========== p <-> q Inversion consistency ===================== > ok <- 1e-5 < pz & pz < 1 - 1e-5 > all.equal(z[ok], qnorm(pz[ok]), tol= 1e-12) [1] TRUE > > ###===== Random numbers -- first, just output: > > .Random.seed <- c(0, 17292, 29447, 24113) > 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.0009459396 0.0749569119 0.7969897349 0.1173196499 0.2400710007 [6] 0.1318767829 0.8277537923 0.0359498142 0.2403329282 0.0372784791 [11] 0.3759311413 0.3581043996 0.2895229973 0.6005864534 0.1833902168 [16] 0.1198140580 0.0706071423 0.0308280998 0.3591918485 0.4281862391 > (Rbinom <- rbinom (n, size = 55, prob = pi/16) ) [1] 12 10 11 13 16 6 5 12 13 9 17 12 14 5 5 13 10 12 16 6 > (Rcauchy <- rcauchy (n, location = 12, scale = 2) ) [1] 34.263018 12.540399 16.802259 11.789490 14.154479 10.462895 12.540308 [8] 12.876666 12.537513 11.977077 11.807015 12.481572 11.954381 10.884467 [15] 9.516357 12.935321 10.946190 14.864709 7.898369 12.782296 > (Rchisq <- rchisq (n, df = 3) ) [1] 4.1802057 2.2410473 8.3668743 0.7165071 4.6012338 0.7650135 2.1202187 [8] 0.2898359 0.7288776 2.9159575 3.6549180 0.5587417 4.9140095 7.3809424 [15] 5.6769453 4.9920360 0.6435048 1.1676611 3.5323618 1.6709330 > (Rexp <- rexp (n, rate = 2) ) [1] 0.07333367 2.13910864 1.32756654 0.11396374 0.07688834 1.38460882 [7] 0.55605730 0.34815444 0.72538091 0.76093410 0.51832937 0.94396435 [13] 0.11381630 0.65756405 0.29282882 0.04916325 0.18753254 0.76833277 [19] 1.01917017 2.40313307 > (Rf <- rf (n, df1 = 12, df2 = 6) ) [1] 1.0482995 7.1506346 0.6029385 2.0461847 1.1851525 4.3366766 0.9365564 [8] 1.2584864 0.4850872 1.4162343 1.3338440 1.1810704 0.7064808 0.5176781 [15] 2.1632007 2.9753123 1.2178294 1.7130207 0.5469131 0.1313706 > (Rgamma <- rgamma (n, shape = 2, scale = 5) ) [1] 3.9092623 21.6048906 11.5498513 11.1312817 4.6122029 14.6707798 [7] 0.9007555 18.9112724 9.4168271 30.4501082 4.5755714 3.2193519 [13] 24.5309039 5.2340905 14.5772699 1.0665363 11.3585753 8.7876520 [19] 8.0414810 14.5654162 > (Rgeom <- rgeom (n, prob = pi/16) ) [1] 4 2 5 1 10 0 0 1 0 17 0 3 1 3 3 6 0 0 2 1 > (Rhyper <- rhyper (n, m = 40, n = 30, k = 20) ) [1] 16 11 11 15 11 13 13 12 13 10 10 7 11 14 13 9 14 13 13 11 > (Rlnorm <- rlnorm (n, meanlog = -1, sdlog = 3) ) [1] 2.210483e-01 1.713584e-01 2.116665e+00 8.205068e+01 3.111700e+00 [6] 1.883586e+00 4.389370e+00 7.965253e-04 3.287647e-01 4.398008e-01 [11] 3.554842e+00 5.186465e-01 2.313136e+00 2.371949e-01 4.293877e-01 [16] 6.577897e-02 2.645768e-01 9.650592e-02 3.650006e-01 7.501445e-02 > (Rlogis <- rlogis (n, location = 12, scale = 2) ) [1] 17.000716 14.220747 22.693644 13.270882 4.748890 16.217577 12.780626 [8] 8.498878 8.300536 9.565866 13.896595 14.425404 17.849150 8.886460 [15] 9.756902 13.182288 15.045885 15.877290 14.773741 9.637214 > (Rnbinom <- rnbinom (n, size = 7, prob = .01) ) [1] 739 668 724 609 1052 404 660 381 321 839 901 292 572 684 883 [16] 1133 503 241 686 787 > (Rnorm <- rnorm (n, mean = -1, sd = 3) ) [1] -0.61029370 2.35074871 2.49862518 0.28586195 1.56824133 4.02175812 [7] -3.40412993 -2.65797756 -0.49471250 -1.78043063 -0.54327159 4.38579193 [13] 0.90049087 -0.06666040 -0.03921922 4.00360536 2.36688992 1.44720283 [19] 1.12529029 0.14489413 > (Rpois <- rpois (n, lambda = 12) ) [1] 11 14 13 7 12 9 19 8 9 12 10 17 11 13 11 12 14 8 6 14 > (Rsignrank<- rsignrank(n, n = 47) ) [1] 521 680 465 368 636 497 560 439 541 462 459 595 542 699 725 563 419 508 774 [20] 496 > (Rt <- rt (n, df = 11) ) [1] -1.73939322 -0.53416778 1.61643746 -0.09307946 0.19731802 0.13339707 [7] -0.99463560 -1.49665151 1.44899521 -0.72605060 -0.52184248 1.62334863 [13] -1.07413952 0.18925757 -0.08940649 -0.90018242 -0.61911928 1.17444416 [19] -1.06421053 -0.77380129 > ## Rt2 below (to preserve the following random numbers!) > (Runif <- runif (n, min = .2, max = 2) ) [1] 1.8416169 0.4703079 0.6841623 0.5048212 1.5620335 0.2738142 1.1573172 [8] 0.3854089 1.0601866 1.9827946 0.7149722 1.3607511 1.3590662 1.3110552 [15] 1.5017535 1.0945284 1.7290782 0.6057555 1.7500334 1.7121475 > (Rweibull <- rweibull (n, shape = 3, scale = 2) ) [1] 1.512345 2.645024 1.458062 2.484227 1.263118 2.081812 1.058986 1.104576 [9] 1.644986 1.985632 2.127461 2.058949 2.524137 1.716619 1.875699 1.433992 [17] 2.362451 2.861329 2.150789 1.239032 > (Rwilcox <- rwilcox (n, m = 13, n = 17) ) [1] 133 103 103 62 117 108 96 102 118 115 110 92 139 126 77 65 127 72 139 [20] 124 > (Rt2 <- rt (n, df = 1.01)) [1] -0.19162182 -0.80558213 5.34371153 2.97056448 -2.06158935 1.02975367 [7] 0.02712426 -1.18826124 1.43878831 2.74129680 -0.52293357 -0.17152843 [13] -0.45179010 1.05646571 -0.62321076 -0.12553755 0.36695132 -9.42934145 [19] 0.04081237 -3.70972593 > > (Pbeta <- pbeta (Rbeta, shape1 = .8, shape2 = 2) ) [1] 0.00685142 0.21898139 0.96943981 0.30726448 0.51350239 0.33510112 [7] 0.97810553 0.12383254 0.51388356 0.12740407 0.68542852 0.66557175 [13] 0.58183279 0.87756637 0.42564683 0.31211434 0.20917242 0.10975854 [19] 0.66680467 0.73943522 > (Pbinom <- pbinom (Rbinom, size = 55, prob = pi/16) ) [1] 0.72571464 0.47335812 0.60646826 0.82208269 0.96870000 0.06545671 [7] 0.02841293 0.72571464 0.82208269 0.34018132 0.98484929 0.72571464 [13] 0.89271720 0.02841293 0.02841293 0.82208269 0.47335812 0.72571464 [19] 0.96870000 0.06545671 > (Pcauchy <- pcauchy (Rcauchy, location = 12, scale = 2) ) [1] 0.9714812 0.5840012 0.8743873 0.4666191 0.7618305 0.2914213 0.5839878 [8] 0.6314968 0.5835731 0.4963519 0.4693803 0.5752128 0.4927407 0.3380482 [15] 0.2157967 0.6392424 0.3456390 0.8059952 0.1444131 0.6186825 > (Pchisq <- pchisq (Rchisq, df = 3) ) [1] 0.75734885 0.47609135 0.96099212 0.13068538 0.79656374 0.14218428 [7] 0.45216677 0.03807139 0.13360847 0.59523480 0.69878142 0.09419094 [13] 0.82179683 0.93930223 0.87156985 0.82761876 0.11359536 0.23922949 [19] 0.68341090 0.35658474 > (Pexp <- pexp (Rexp, rate = 2) ) [1] 0.13641880 0.98613264 0.92971052 0.20381800 0.14253652 0.93728895 [7] 0.67113719 0.50157835 0.76560834 0.78169633 0.64536236 0.84861495 [13] 0.20358320 0.73156007 0.44326040 0.09364707 0.31275544 0.78490286 [19] 0.86975531 0.99182166 > (Pf <- pf (Rf, df1 = 12, df2 = 6) ) [1] 0.493968892 0.987819917 0.214465434 0.804851021 0.560142915 0.958476676 [7] 0.432468011 0.591675896 0.134821391 0.651135448 0.621437022 0.558309651 [13] 0.285976592 0.156229064 0.823256086 0.904828024 0.574518093 0.737314198 [19] 0.175931308 0.001535476 > (Pgamma <- pgamma (Rgamma, shape = 2, scale = 5) ) [1] 0.18469983 0.92930080 0.67143844 0.65178456 0.23573834 0.79079896 [7] 0.01440388 0.89110187 0.56150517 0.98393855 0.23305255 0.13654797 [13] 0.95629013 0.28145838 0.78786294 0.01975953 0.66257119 0.52440798 [19] 0.47774483 0.78748817 > (Pgeom <- pgeom (Rgeom, prob = pi/16) ) [1] 0.6647753 0.4809591 0.7305965 0.3541459 0.9096893 0.1963495 0.1963495 [8] 0.3541459 0.1963495 0.9804472 0.1963495 0.5828725 0.3541459 0.5828725 [15] 0.5828725 0.7834938 0.1963495 0.1963495 0.4809591 0.3541459 > (Phyper <- phyper (Rhyper, m = 40, n = 30, k = 20) ) [1] 0.99744478 0.51295659 0.51295659 0.98680472 0.51295659 0.86627413 [7] 0.86627413 0.71494883 0.86627413 0.30864260 0.30864260 0.01796063 [13] 0.51295659 0.95139461 0.86627413 0.15132082 0.95139461 0.86627413 [19] 0.86627413 0.51295659 > (Plnorm <- plnorm (Rlnorm, meanlog = -1, sdlog = 3) ) [1] 0.43258714 0.39949055 0.72014780 0.96426235 0.76168186 0.70691429 [7] 0.79571003 0.02042329 0.48505472 0.52373190 0.77520644 0.54557496 [13] 0.73001823 0.44184582 0.52055037 0.28304536 0.45625441 0.32778045 [19] 0.49895528 0.29804735 > (Plogis <- plogis (Rlogis, location = 12, scale = 2) ) [1] 0.92416690 0.75219874 0.99525932 0.65372216 0.02594333 0.89175446 [7] 0.59635809 0.14797648 0.13590439 0.22845301 0.72077271 0.77077670 [13] 0.94904799 0.17411062 0.24572404 0.64362757 0.82097136 0.87420320 [19] 0.80009214 0.23480185 > (Pnbinom <- pnbinom (Rnbinom, size = 7, prob = .01) ) [1] 0.61736655 0.51299455 0.59632006 0.41950678 0.90402573 0.12140233 [7] 0.50057859 0.09750529 0.04876502 0.74086686 0.80178090 0.03225449 [13] 0.35983768 0.53748456 0.78532898 0.93731909 0.25190327 0.01316344 [19] 0.54051099 0.68045166 > (Pnorm <- pnorm (Rnorm, mean = -1, sd = 3) ) [1] 0.5516781 0.8679849 0.8782349 0.6659003 0.8040234 0.9529268 0.2114568 [8] 0.2902484 0.5668772 0.3973765 0.5605023 0.9636942 0.7367954 0.6221427 [15] 0.6256145 0.9523291 0.8691318 0.7926739 0.7606609 0.6486326 > (Ppois <- ppois (Rpois, lambda = 12) ) [1] 0.46159733 0.77202453 0.68153563 0.08950450 0.57596525 0.24239216 [7] 0.97872023 0.15502778 0.24239216 0.57596525 0.34722942 0.93703370 [13] 0.46159733 0.68153563 0.46159733 0.57596525 0.77202453 0.15502778 [19] 0.04582231 0.77202453 > (Psignrank<- psignrank(Rsignrank, n = 47) ) [1] 0.32789203 0.89007439 0.15011688 0.01887176 0.77676020 0.24256054 [7] 0.48536796 0.09474662 0.40676506 0.14282568 0.13577317 0.62940572 [13] 0.41083885 0.92356606 0.95623357 0.49790926 0.06352849 0.28016760 [19] 0.98752692 0.23928001 > (Pt <- pt (Rt, df = 11) ) [1] 0.05491642 0.30192660 0.93285384 0.46375705 0.57641129 0.55185574 [7] 0.17064488 0.08130751 0.91237971 0.24148527 0.30606307 0.93360050 [13] 0.15287884 0.57333111 0.46518298 0.19365307 0.27422273 0.86749203 [19] 0.15501842 0.22767857 > (Pt2 <- pt (Rt2, df = 1.01) ) [1] 0.43961636 0.28360180 0.94199308 0.89763502 0.14278894 0.75533914 [7] 0.50864836 0.22193161 0.80750963 0.88965697 0.34628392 0.44582080 [13] 0.36463226 0.75942529 0.32218488 0.46017088 0.61218435 0.03294746 [19] 0.51300870 0.08284141 > (Punif <- punif (Runif, min = .2, max = 2) ) [1] 0.91200940 0.15017105 0.26897908 0.16934513 0.75668526 0.04100786 [7] 0.53184287 0.10300496 0.47788145 0.99044142 0.28609568 0.64486175 [13] 0.64392568 0.61725289 0.72319637 0.49696022 0.84948787 0.22541974 [19] 0.86112968 0.84008195 > (Pweibull <- pweibull (Rweibull, shape = 3, scale = 2) ) [1] 0.3510356 0.9010483 0.3212279 0.8528629 0.2226838 0.6762575 0.1379570 [8] 0.1550349 0.4267381 0.6241639 0.6998987 0.6641379 0.8660439 0.4686384 [15] 0.5617197 0.3082951 0.8075955 0.9465114 0.7116742 0.2116163 > (Pwilcox <- pwilcox (Rwilcox, m = 13, n = 17) ) [1] 0.82946462 0.38647276 0.38647276 0.02165645 0.61352724 0.46716310 [7] 0.28170171 0.37078117 0.62921883 0.58162230 0.50000000 0.22852300 [13] 0.88566977 0.74559308 0.08487632 0.02943715 0.75871780 0.05633406 [19] 0.88566977 0.71829829 > > dbeta (Rbeta, shape1 = .8, shape2 = 2) [1] 5.7913365 2.2364651 0.3059074 1.9511551 1.4556855 1.8746068 0.2575917 [8] 2.6997454 1.4548665 2.6765264 1.0928822 1.1350760 1.3109191 0.6368984 [15] 1.6508280 1.9374717 2.2740086 2.7988135 1.1324661 0.9756421 > dbinom (Rbinom, size = 55, prob = pi/16) [1] 0.11924639 0.13317680 0.13311013 0.09636804 0.02881210 0.03704378 [7] 0.01819424 0.11924639 0.09636804 0.11849719 0.01614929 0.11924639 [13] 0.07063452 0.01819424 0.01819424 0.09636804 0.13317680 0.11924639 [19] 0.02881210 0.03704378 > dcauchy (Rcauchy, location = 12, scale = 2) [1] 0.001274152 0.148326009 0.023524757 0.157411041 0.073667654 0.100055114 [7] 0.148329389 0.133504054 0.148433568 0.159134039 0.157686757 0.150433150 [13] 0.159072180 0.121390049 0.062607155 0.130593337 0.124570555 0.052153921 [19] 0.030572357 0.138035910 > dchisq (Rchisq, df = 3) [1] 0.10087621 0.19475936 0.01759333 0.23601129 0.08574365 0.23802581 [7] 0.20123372 0.18580144 0.23657215 0.15852893 0.12265734 0.22552069 [13] 0.07578172 0.02705295 0.05562044 0.07345850 0.23198007 0.24044408 [19] 0.12820355 0.22364083 > dexp (Rexp, rate = 2) [1] 1.72716241 0.02773472 0.14057897 1.59236400 1.71492696 0.12542210 [7] 0.65772563 0.99684330 0.46878331 0.43660735 0.70927528 0.30277010 [13] 1.59283360 0.53687987 1.11347921 1.81270585 1.37448911 0.43019428 [19] 0.26048939 0.01635668 > df (Rf, df1 = 12, df2 = 6) [1] 0.519937867 0.004371848 0.692823649 0.167468952 0.448061444 0.022237248 [7] 0.580744146 0.412297607 0.645723076 0.343445106 0.377992320 0.450116704 [13] 0.682397435 0.666844216 0.147593859 0.066212131 0.431848359 0.243311742 [19] 0.680166477 0.051540741 > dgamma (Rgamma, shape = 2, scale = 5) [1] 0.071548510 0.011482464 0.045859471 0.048056741 0.073342809 0.031205091 [7] 0.030090422 0.017225367 0.057283478 0.002759218 0.073295320 0.067639298 [13] 0.007261834 0.073497725 0.031591526 0.034466522 0.046858733 0.060624349 [19] 0.064405338 0.031640760 > dgeom (Rgeom, prob = pi/16) [1] 0.081902787 0.126813148 0.065821212 0.157796399 0.022064895 0.196349541 [7] 0.196349541 0.157796399 0.196349541 0.004777187 0.196349541 0.101913445 [13] 0.157796399 0.101913445 0.101913445 0.052897247 0.196349541 0.196349541 [19] 0.126813148 0.157796399 > dhyper (Rhyper, m = 40, n = 30, k = 20) [1] 0.01064006 0.20431399 0.20431399 0.03541012 0.20431399 0.15132529 [7] 0.15132529 0.20199224 0.15132529 0.15732178 0.15732178 0.01379223 [13] 0.20431399 0.08512048 0.15132529 0.09227084 0.08512048 0.15132529 [19] 0.15132529 0.20431399 > dlnorm (Rlnorm, meanlog = -1, sdlog = 3) [1] 5.929818e-01 7.512774e-01 5.299802e-02 3.193275e-04 3.317380e-02 [6] 6.087639e-02 2.153230e-02 2.062527e+01 4.042023e-01 3.018308e-01 [11] 2.810781e-02 2.547247e-01 4.764589e-02 5.546720e-01 3.092876e-01 [16] 1.714756e+00 4.995922e-01 1.247473e+00 3.643290e-01 1.540428e+00 > dlogis (Rlogis, location = 12, scale = 2) [1] 0.035041221 0.093197897 0.002359104 0.113184749 0.012635138 0.048264223 [7] 0.120357559 0.063039721 0.058717192 0.088131115 0.100629705 0.088339990 [13] 0.024177951 0.071898055 0.092671867 0.114685561 0.073488693 0.054985983 [19] 0.079972353 0.089834971 > dnbinom (Rnbinom, size = 7, prob = .01) [1] 0.0013843572 0.0015461270 0.0014240958 0.0016111183 0.0004916298 [6] 0.0010965512 0.0015593115 0.0009750548 0.0006438980 0.0010814585 [11] 0.0008880238 0.0004914120 0.0016079010 0.0015162429 0.0009432177 [16] 0.0003394350 0.0014949693 0.0002632443 0.0015121979 0.0012444310 > dnorm (Rnorm, mean = -1, sd = 3) [1] 0.13186348 0.07126843 0.06736890 0.12130968 0.09218225 0.03275978 [7] 0.09645747 0.11414715 0.13110785 0.12855632 0.13144855 0.02654173 [13] 0.10880369 0.12669831 0.12633298 0.03309267 0.07084041 0.09534418 [19] 0.10346857 0.12364112 > dpois (Rpois, lambda = 12) [1] 0.11436792 0.09048890 0.10557038 0.04368219 0.11436792 0.08736438 [7] 0.01613672 0.06552328 0.08736438 0.11436792 0.10483726 0.03832471 [13] 0.11436792 0.10557038 0.11436792 0.11436792 0.09048890 0.06552328 [19] 0.02548128 0.09048890 > dsignrank(Rsignrank, n = 47) [1] 0.0037844241 0.0020112978 0.0024569596 0.0004983669 0.0031591482 [6] 0.0032805288 0.0041778794 0.0017851857 0.0040639142 0.0023773231 [11] 0.0022979400 0.0039702993 0.0040737886 0.0015466398 0.0010084031 [16] 0.0041812581 0.0013244780 0.0035300054 0.0003600738 0.0032566009 > dt (Rt, df = 11) [1] 0.09076142 0.33444332 0.10857032 0.38815185 0.38180921 0.38622576 [7] 0.23261963 0.12825909 0.13672995 0.29449489 0.33676702 0.10750499 [13] 0.21435914 0.38245649 0.38829368 0.25458773 0.31753985 0.19196797 [19] 0.21661975 0.28375142 > 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.5566140 0.2596054 0.5411369 0.3405149 0.4650686 0.5261555 0.3625272 [8] 0.3865994 0.5817133 0.5556831 0.5093571 0.5339287 0.3200512 0.5871773 [15] 0.5782419 0.5333911 0.4026909 0.1642208 0.5001614 0.4538727 > dwilcox (Rwilcox, m = 13, n = 17) [1] 0.010799421 0.015691594 0.015691594 0.002185415 0.015873676 0.016351073 [7] 0.013817085 0.015486075 0.015691594 0.016163798 0.016431951 0.012383032 [13] 0.008348499 0.013477697 0.006413577 0.002809857 0.013124724 0.004704515 [19] 0.008348499 0.014141025 > > ## 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.equal(Rt2, qt (Pt2, df = 1.01), tol = 1e-2) [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.equal(Rt2, qt (1- Pt2, df = 1.01, lower=F), tol = 1e-2) [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.equal(Rchisq, qchisq (log(Pchisq), df = 3, log=TRUE),tol=1e-14) [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.equal(Rt2, qt (log(Pt2), df = 1.01, log=TRUE), tol = 1e-2) [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 (log(1- Pbeta), shape1 = .8, shape2 = 2, lower=F, log=T)) [1] TRUE > All.eq(Rbinom, qbinom (log(1- Pbinom), size = 55, prob = pi/16, lower=F, log=T)) [1] TRUE > All.eq(Rcauchy, qcauchy (log(1- Pcauchy), location = 12, scale = 2, lower=F, log=T)) [1] TRUE > All.eq(Rchisq, qchisq (log(1- Pchisq), df = 3, lower=F, log=T)) [1] TRUE > All.eq(Rexp, qexp (log(1- Pexp), rate = 2, lower=F, log=T)) [1] TRUE > All.eq(Rf, qf (log(1- Pf), df1 = 12, df2 = 6, lower=F, log=T)) [1] TRUE > All.eq(Rgamma, qgamma (log(1- Pgamma), shape = 2, scale = 5, lower=F, log=T)) [1] TRUE > All.eq(Rgeom, qgeom (log(1- Pgeom), prob = pi/16, lower=F, log=T)) [1] TRUE > All.eq(Rhyper, qhyper (log(1- Phyper), m = 40, n = 30, k = 20, lower=F, log=T)) [1] TRUE > All.eq(Rlnorm, qlnorm (log(1- Plnorm), meanlog = -1, sdlog = 3, lower=F, log=T)) [1] TRUE > All.eq(Rlogis, qlogis (log(1- Plogis), location = 12, scale = 2, lower=F, log=T)) [1] TRUE > All.eq(Rnbinom, qnbinom (log(1- Pnbinom), size = 7, prob = .01, lower=F, log=T)) [1] TRUE > All.eq(Rnorm, qnorm (log(1- Pnorm), mean = -1, sd = 3, lower=F, log=T)) [1] TRUE > All.eq(Rpois, qpois (log(1- Ppois), lambda = 12, lower=F, log=T)) [1] TRUE > All.eq(Rsignrank, qsignrank(log(1- Psignrank), n = 47, lower=F, log=T)) [1] TRUE > All.eq(Rt, qt (log(1- Pt ), df = 11, lower=F, log=T)) [1] TRUE > all.equal(Rt2, qt (log(1- Pt2), df = 1.01, lower=F, log=T), tol = 1e-2) [1] TRUE > All.eq(Runif, qunif (log(1- Punif), min = .2, max = 2, lower=F, log=T)) [1] TRUE > All.eq(Rweibull, qweibull (log(1- Pweibull), shape = 3, scale = 2, lower=F, log=T)) [1] TRUE > All.eq(Rwilcox, qwilcox (log(1- Pwilcox), m = 13, n = 17, lower=F, log=T)) [1] TRUE > > > ## Check log( upper.tail ): > All.eq(log(1 - Pbeta), pbeta (Rbeta, shape1 = .8, shape2 = 2, lower=F, log=T)) [1] TRUE > All.eq(log(1 - Pbinom), pbinom (Rbinom, size = 55, prob = pi/16, lower=F, log=T)) [1] TRUE > All.eq(log(1 - Pcauchy), pcauchy (Rcauchy, location = 12, scale = 2, lower=F, log=T)) [1] TRUE > All.eq(log(1 - Pchisq), pchisq (Rchisq, df = 3, lower=F, log=T)) [1] TRUE > All.eq(log(1 - Pexp), pexp (Rexp, rate = 2, lower=F, log=T)) [1] TRUE > All.eq(log(1 - Pf), pf (Rf, df1 = 12, df2 = 6, lower=F, log=T)) [1] TRUE > All.eq(log(1 - Pgamma), pgamma (Rgamma, shape = 2, scale = 5, lower=F, log=T)) [1] TRUE > All.eq(log(1 - Pgeom), pgeom (Rgeom, prob = pi/16, lower=F, log=T)) [1] TRUE > All.eq(log(1 - Phyper), phyper (Rhyper, m = 40, n = 30, k = 20, lower=F, log=T)) [1] TRUE > All.eq(log(1 - Plnorm), plnorm (Rlnorm, meanlog = -1, sdlog = 3, lower=F, log=T)) [1] TRUE > All.eq(log(1 - Plogis), plogis (Rlogis, location = 12, scale = 2, lower=F, log=T)) [1] TRUE > All.eq(log(1 - Pnbinom), pnbinom (Rnbinom, size = 7, prob = .01, lower=F, log=T)) [1] TRUE > All.eq(log(1 - Pnorm), pnorm (Rnorm, mean = -1, sd = 3, lower=F, log=T)) [1] TRUE > All.eq(log(1 - Ppois), ppois (Rpois, lambda = 12, lower=F, log=T)) [1] TRUE > All.eq(log(1 - Psignrank), psignrank(Rsignrank, n = 47, lower=F, log=T)) [1] TRUE > All.eq(log(1 - Pt), pt (Rt, df = 11, lower=F, log=T)) [1] TRUE > All.eq(log(1 - Pt2), pt (Rt2,df = 1.01, lower=F, log=T)) [1] TRUE > All.eq(log(1 - Punif), punif (Runif, min = .2, max = 2, lower=F, log=T)) [1] TRUE > All.eq(log(1 - Pweibull), pweibull (Rweibull, shape = 3, scale = 2, lower=F, log=T)) [1] TRUE > All.eq(log(1 - 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)) + print(all.equal(pt(-a, df=1), pcauchy(-a), tol = 1e-15)) [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE [1] TRUE > ## 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 <- 1e-23 ## PR#6757 > stopifnot(all.equal(pr^ 12, pbinom(11, 12, prob= pr,lower=FALSE), + tol= 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)), tol = 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)), tol = 1e-13)# 7e-16 (i386) [1] TRUE > > All.eq(dpois( 10, 2e-308, log=TRUE), -7100.13502718914) [1] TRUE > All.eq(dpois( 20, 3e-308, log=TRUE), -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 > pf(x, 5, 1e6, ncp=1) [1] 0.065933194 0.470879987 0.978875867 > pf(x, 5, 1e7, ncp=1) [1] 0.06593309 0.47088028 0.97887641 > all.equal(pf(x, 5, 1e8, ncp=1), tol = 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) > > ## PR#7099 : pf() with large df1 or df2: > nu <- 2^seq(25,34, 0.5) > y <- 1e9*(pf(1,1,nu) - 0.68268949) > stopifnot(All.eq(pf(1,1,Inf), 0.68268949213708596), + diff(y) > 0, # i.e. pf(1,1, *) is monotone increasing + All.eq(y [1], -5.07420372386491), + All.eq(y[19], 2.12300110824515)) > ## not at all in R 2.1.0 or earlier > > 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 > > 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) > stopifnot(pbinom(x0, size = 3, prob = 0.1) == 0, + dbinom(x0, 3, 0.1) == 0) # d*() warns about non-integer Warning messages: 1: non-integer x = -0.000000 2: non-integer x = -0.000020 > ## very small negatives were rounded to 0 in R 2.2.1 and earlier > > ## dbeta(*, ncp): > a <- rlnorm(100) > stopifnot(All.eq(a, dbeta(0, 1, a, ncp=0)), + dbeta(0, 0.9, 2.2, ncp = c(0, a)) == Inf + ) > ## the first gave 0, the 2nd NaN in R <= 2.3.0 > > ## 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), tol = 1e-6),# (1.50 | 1.65)e-8 + All.eq(df(0, 2, 4, ncp=x), df(1e-300, 2, 4, ncp=x)) + ) > > cat("Time elapsed: ", proc.time() - .ptime,"\n") Time elapsed: 5.148 0.034 5.203 0 0 >