R : Copyright 2001, The R Development Core Team
Version 1.4.0 Under development (unstable) (2001-09-07)
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.
Type `demo()' for some demos, `help()' for on-line help, or
`help.start()' for a 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)
+ }
> ## Short cut:
> All.eq <- function(x,y) all.equal.numeric(x,y, tolerance = 1e-14)
>
> 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(!(is.logical(tst1) && tst1) ||
+ !(is.logical(tst) && tst ) ) {
+ cat("n=", n,"; p =",format(p),". k =",k)
+ if(!is.logical(tst1)) cat("; tst1=",tst1)
+ if(!is.logical(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(0,1,len=15)) {
+ 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(!(is.logical(tst) && tst))
+ cat("lambda=", format(lambda),". k =",k, " --> tst=", tst,"\n")
+ tst2 <- all.equal(pp, ppois(0:k, lambda=lambda), tol = 100*Meps)
+ if(!(is.logical(tst2) && tst2))
+ cat("lambda=", format(lambda),". k =",k, " --> tst2=", tst2,"\n")
+ tst3 <- all.equal(1 - pp, ppois(0:k, lambda=lambda, lower.tail=FALSE))
+ if(!(is.logical(tst3) && tst3))
+ cat("lambda=", format(lambda),". k =",k, " --> tst3=", tst3,"\n")
+ }
>
>
> ##__ 6. SignRank __
> for(n in rpois(32, lam=8)) {
+ x <- -1:(n + 4)
+ eq <- All.eq(psignrank(x, n), cumsum(dsignrank(x, n)))
+ if(!is.logical(eq) || !eq) 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))
+ eq <- All.eq(Fx, cumsum(fx))
+ if(!is.logical(eq) || !eq) print(eq)
+ }
> is.sym
[1] TRUE
>
>
> ###-------- Continuous Distributions ----------
>
> ##--- Gamma (incl. 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(!(is.logical(tst) && 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(!(is.logical(tst) && 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
> all(is.nan(c(pgamma(Inf,1,scale=Inf), pgamma(Inf,Inf,scale=1), 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)
3: NaNs produced in: pgamma(q, shape, scale, lower.tail, log.p)
> pgamma(Inf,1,scale=xMax) == 1 && pgamma(xMax,1,scale=Inf) == 0
[1] TRUE
>
> ##--- 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)
>
> ## 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(!is.logical(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
>
> 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
> (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
>
> (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
> (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.equal( Rt, qt (Pt, df = 11), tol = 1e-7 )## !! (386-lx:2.48e-9)
[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.equal(Rt, qt (1- Pt, df = 11, lower=F), tol = 1e-7)
[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.equal(Rt, qt (log(Pt), df = 11, log=TRUE), tol = 1e-7)
[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.equal(Rt, qt (log(1- Pt), df = 11, lower=F, log=T), tol = 1e-7)#ok{prec.}
[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 - 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
>
> cat("Time elapsed: ", proc.time() - .ptime,"\n")
Time elapsed: 44.99 7.5 53.2 0 0
>