#### 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 = 100*.Machine$double.eps) 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") } ##__ 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))) } ##__ 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)))) } ##__ 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))) All.eq(pnbinom(c(1,3), .9, .5), c(0.777035760338812, 0.946945347071519)) ##__ 5. Poisson __ all(dpois(0:5,0) == c(1, rep(0,5))) all(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)) { 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 ###-------- 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(!(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 all(is.nan(c(pgamma(Inf,1,scale=Inf), pgamma(Inf,Inf,scale=1), pgamma(Inf,Inf,scale=Inf)))) pgamma(Inf,1,scale=xMax) == 1 && pgamma(xMax,1,scale=Inf) == 0 ##-- 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) ##--- Beta (need more): ## big a & b (PR #643) summary(a <- rlnorm(20, 5.5)) summary(b <- rlnorm(20, 6.5)) 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) sample(pab, 50) ##--- Normal (& Lognormal) : qnorm(0) == -Inf && qnorm(-Inf, log = TRUE) == -Inf qnorm(1) == Inf && qnorm(0, log = TRUE) == Inf is.nan(qnorm(1.1)) && is.nan(qnorm(-.1)) # + warn ## 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) # 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), tol= 1e-15) 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)) 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 <- 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)) ## 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))) ###========== p <-> q Inversion consistency ===================== ok <- 1e-5 < pz & pz < 1 - 1e-5 all.equal(z[ok], qnorm(pz[ok]), tol= 1e-12) ###===== 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) ) (Rbinom <- rbinom (n, size = 55, prob = pi/16) ) (Rcauchy <- rcauchy (n, location = 12, scale = 2) ) (Rchisq <- rchisq (n, df = 3) ) (Rexp <- rexp (n, rate = 2) ) (Rf <- rf (n, df1 = 12, df2 = 6) ) (Rgamma <- rgamma (n, shape = 2, scale = 5) ) (Rgeom <- rgeom (n, prob = pi/16) ) (Rhyper <- rhyper (n, m = 40, n = 30, k = 20) ) (Rlnorm <- rlnorm (n, meanlog = -1, sdlog = 3) ) (Rlogis <- rlogis (n, location = 12, scale = 2) ) (Rnbinom <- rnbinom (n, size = 7, prob = .01) ) (Rnorm <- rnorm (n, mean = -1, sd = 3) ) (Rpois <- rpois (n, lambda = 12) ) (Rsignrank<- rsignrank(n, n = 47) ) (Rt <- rt (n, df = 11) ) ## Rt2 below (to preserve the following random numbers!) (Runif <- runif (n, min = .2, max = 2) ) (Rweibull <- rweibull (n, shape = 3, scale = 2) ) (Rwilcox <- rwilcox (n, m = 13, n = 17) ) (Rt2 <- rt (n, df = 1.01)) (Pbeta <- pbeta (Rbeta, shape1 = .8, shape2 = 2) ) (Pbinom <- pbinom (Rbinom, size = 55, prob = pi/16) ) (Pcauchy <- pcauchy (Rcauchy, location = 12, scale = 2) ) (Pchisq <- pchisq (Rchisq, df = 3) ) (Pexp <- pexp (Rexp, rate = 2) ) (Pf <- pf (Rf, df1 = 12, df2 = 6) ) (Pgamma <- pgamma (Rgamma, shape = 2, scale = 5) ) (Pgeom <- pgeom (Rgeom, prob = pi/16) ) (Phyper <- phyper (Rhyper, m = 40, n = 30, k = 20) ) (Plnorm <- plnorm (Rlnorm, meanlog = -1, sdlog = 3) ) (Plogis <- plogis (Rlogis, location = 12, scale = 2) ) (Pnbinom <- pnbinom (Rnbinom, size = 7, prob = .01) ) (Pnorm <- pnorm (Rnorm, mean = -1, sd = 3) ) (Ppois <- ppois (Rpois, lambda = 12) ) (Psignrank<- psignrank(Rsignrank, n = 47) ) (Pt <- pt (Rt, df = 11) ) (Pt2 <- pt (Rt2, df = 1.01) ) (Punif <- punif (Runif, min = .2, max = 2) ) (Pweibull <- pweibull (Rweibull, shape = 3, scale = 2) ) (Pwilcox <- pwilcox (Rwilcox, m = 13, n = 17) ) dbeta (Rbeta, shape1 = .8, shape2 = 2) dbinom (Rbinom, size = 55, prob = pi/16) dcauchy (Rcauchy, location = 12, scale = 2) dchisq (Rchisq, df = 3) dexp (Rexp, rate = 2) df (Rf, df1 = 12, df2 = 6) dgamma (Rgamma, shape = 2, scale = 5) dgeom (Rgeom, prob = pi/16) dhyper (Rhyper, m = 40, n = 30, k = 20) dlnorm (Rlnorm, meanlog = -1, sdlog = 3) dlogis (Rlogis, location = 12, scale = 2) dnbinom (Rnbinom, size = 7, prob = .01) dnorm (Rnorm, mean = -1, sd = 3) dpois (Rpois, lambda = 12) dsignrank(Rsignrank, n = 47) dt (Rt, df = 11) dunif (Runif, min = .2, max = 2) dweibull (Rweibull, shape = 3, scale = 2) dwilcox (Rwilcox, m = 13, n = 17) ## Check q*(p*(.)) = identity All.eq(Rbeta, qbeta (Pbeta, shape1 = .8, shape2 = 2)) All.eq(Rbinom, qbinom (Pbinom, size = 55, prob = pi/16)) All.eq(Rcauchy, qcauchy (Pcauchy, location = 12, scale = 2)) All.eq(Rchisq, qchisq (Pchisq, df = 3)) All.eq(Rexp, qexp (Pexp, rate = 2)) All.eq(Rf, qf (Pf, df1 = 12, df2 = 6)) All.eq(Rgamma, qgamma (Pgamma, shape = 2, scale = 5)) All.eq(Rgeom, qgeom (Pgeom, prob = pi/16)) All.eq(Rhyper, qhyper (Phyper, m = 40, n = 30, k = 20)) All.eq(Rlnorm, qlnorm (Plnorm, meanlog = -1, sdlog = 3)) All.eq(Rlogis, qlogis (Plogis, location = 12, scale = 2)) All.eq(Rnbinom, qnbinom (Pnbinom, size = 7, prob = .01)) All.eq(Rnorm, qnorm (Pnorm, mean = -1, sd = 3)) All.eq(Rpois, qpois (Ppois, lambda = 12)) All.eq(Rsignrank, qsignrank(Psignrank, n = 47)) All.eq(Rt, qt (Pt, df = 11)) all.equal(Rt2, qt (Pt2, df = 1.01), tol = 1e-2) All.eq(Runif, qunif (Punif, min = .2, max = 2)) All.eq(Rweibull, qweibull (Pweibull, shape = 3, scale = 2)) All.eq(Rwilcox, qwilcox (Pwilcox, m = 13, n = 17)) ## Same with "upper tail": All.eq(Rbeta, qbeta (1- Pbeta, shape1 = .8, shape2 = 2, lower=F)) All.eq(Rbinom, qbinom (1- Pbinom, size = 55, prob = pi/16, lower=F)) All.eq(Rcauchy, qcauchy (1- Pcauchy, location = 12, scale = 2, lower=F)) All.eq(Rchisq, qchisq (1- Pchisq, df = 3, lower=F)) All.eq(Rexp, qexp (1- Pexp, rate = 2, lower=F)) All.eq(Rf, qf (1- Pf, df1 = 12, df2 = 6, lower=F)) All.eq(Rgamma, qgamma (1- Pgamma, shape = 2, scale = 5, lower=F)) All.eq(Rgeom, qgeom (1- Pgeom, prob = pi/16, lower=F)) All.eq(Rhyper, qhyper (1- Phyper, m = 40, n = 30, k = 20, lower=F)) All.eq(Rlnorm, qlnorm (1- Plnorm, meanlog = -1, sdlog = 3, lower=F)) All.eq(Rlogis, qlogis (1- Plogis, location = 12, scale = 2, lower=F)) All.eq(Rnbinom, qnbinom (1- Pnbinom, size = 7, prob = .01, lower=F)) All.eq(Rnorm, qnorm (1- Pnorm, mean = -1, sd = 3,lower=F)) All.eq(Rpois, qpois (1- Ppois, lambda = 12, lower=F)) All.eq(Rsignrank, qsignrank(1- Psignrank, n = 47, lower=F)) All.eq(Rt, qt (1- Pt, df = 11, lower=F)) all.equal(Rt2, qt (1- Pt2, df = 1.01, lower=F), tol = 1e-2) All.eq(Runif, qunif (1- Punif, min = .2, max = 2, lower=F)) All.eq(Rweibull, qweibull (1- Pweibull, shape = 3, scale = 2, lower=F)) All.eq(Rwilcox, qwilcox (1- Pwilcox, m = 13, n = 17, lower=F)) ## Check q*(p* ( log ), log) = identity All.eq(Rbeta, qbeta (log(Pbeta), shape1 = .8, shape2 = 2, log=TRUE)) All.eq(Rbinom, qbinom (log(Pbinom), size = 55, prob = pi/16, log=TRUE)) All.eq(Rcauchy, qcauchy (log(Pcauchy), location = 12, scale = 2, log=TRUE)) all.equal(Rchisq, qchisq (log(Pchisq), df = 3, log=TRUE),tol=1e-14) All.eq(Rexp, qexp (log(Pexp), rate = 2, log=TRUE)) All.eq(Rf, qf (log(Pf), df1= 12, df2= 6, log=TRUE)) All.eq(Rgamma, qgamma (log(Pgamma), shape = 2, scale = 5, log=TRUE)) All.eq(Rgeom, qgeom (log(Pgeom), prob = pi/16, log=TRUE)) All.eq(Rhyper, qhyper (log(Phyper), m = 40, n = 30, k = 20, log=TRUE)) All.eq(Rlnorm, qlnorm (log(Plnorm), meanlog = -1, sdlog = 3, log=TRUE)) All.eq(Rlogis, qlogis (log(Plogis), location = 12, scale = 2, log=TRUE)) All.eq(Rnbinom, qnbinom (log(Pnbinom), size = 7, prob = .01, log=TRUE)) All.eq(Rnorm, qnorm (log(Pnorm), mean = -1, sd = 3, log=TRUE)) All.eq(Rpois, qpois (log(Ppois), lambda = 12, log=TRUE)) All.eq(Rsignrank, qsignrank(log(Psignrank), n = 47, log=TRUE)) All.eq(Rt, qt (log(Pt), df = 11, log=TRUE)) all.equal(Rt2, qt (log(Pt2), df = 1.01, log=TRUE), tol = 1e-2) All.eq(Runif, qunif (log(Punif), min = .2, max = 2, log=TRUE)) All.eq(Rweibull, qweibull (log(Pweibull), shape = 3, scale = 2, log=TRUE)) All.eq(Rwilcox, qwilcox (log(Pwilcox), m = 13, n = 17, log=TRUE)) ## same q*(p* (log) log) with upper tail: All.eq(Rbeta, qbeta (log(1- Pbeta), shape1 = .8, shape2 = 2, lower=F, log=T)) All.eq(Rbinom, qbinom (log(1- Pbinom), size = 55, prob = pi/16, lower=F, log=T)) All.eq(Rcauchy, qcauchy (log(1- Pcauchy), location = 12, scale = 2, lower=F, log=T)) All.eq(Rchisq, qchisq (log(1- Pchisq), df = 3, lower=F, log=T)) All.eq(Rexp, qexp (log(1- Pexp), rate = 2, lower=F, log=T)) All.eq(Rf, qf (log(1- Pf), df1 = 12, df2 = 6, lower=F, log=T)) All.eq(Rgamma, qgamma (log(1- Pgamma), shape = 2, scale = 5, lower=F, log=T)) All.eq(Rgeom, qgeom (log(1- Pgeom), prob = pi/16, lower=F, log=T)) All.eq(Rhyper, qhyper (log(1- Phyper), m = 40, n = 30, k = 20, lower=F, log=T)) All.eq(Rlnorm, qlnorm (log(1- Plnorm), meanlog = -1, sdlog = 3, lower=F, log=T)) All.eq(Rlogis, qlogis (log(1- Plogis), location = 12, scale = 2, lower=F, log=T)) All.eq(Rnbinom, qnbinom (log(1- Pnbinom), size = 7, prob = .01, lower=F, log=T)) All.eq(Rnorm, qnorm (log(1- Pnorm), mean = -1, sd = 3, lower=F, log=T)) All.eq(Rpois, qpois (log(1- Ppois), lambda = 12, lower=F, log=T)) All.eq(Rsignrank, qsignrank(log(1- Psignrank), n = 47, lower=F, log=T)) All.eq(Rt, qt (log(1- Pt ), df = 11, lower=F, log=T)) all.equal(Rt2, qt (log(1- Pt2), df = 1.01, lower=F, log=T), tol = 1e-2) All.eq(Runif, qunif (log(1- Punif), min = .2, max = 2, lower=F, log=T)) All.eq(Rweibull, qweibull (log(1- Pweibull), shape = 3, scale = 2, lower=F, log=T)) All.eq(Rwilcox, qwilcox (log(1- Pwilcox), m = 13, n = 17, lower=F, log=T)) ## Check log( upper.tail ): All.eq(log(1 - Pbeta), pbeta (Rbeta, shape1 = .8, shape2 = 2, lower=F, log=T)) All.eq(log(1 - Pbinom), pbinom (Rbinom, size = 55, prob = pi/16, lower=F, log=T)) All.eq(log(1 - Pcauchy), pcauchy (Rcauchy, location = 12, scale = 2, lower=F, log=T)) All.eq(log(1 - Pchisq), pchisq (Rchisq, df = 3, lower=F, log=T)) All.eq(log(1 - Pexp), pexp (Rexp, rate = 2, lower=F, log=T)) All.eq(log(1 - Pf), pf (Rf, df1 = 12, df2 = 6, lower=F, log=T)) All.eq(log(1 - Pgamma), pgamma (Rgamma, shape = 2, scale = 5, lower=F, log=T)) All.eq(log(1 - Pgeom), pgeom (Rgeom, prob = pi/16, lower=F, log=T)) All.eq(log(1 - Phyper), phyper (Rhyper, m = 40, n = 30, k = 20, lower=F, log=T)) All.eq(log(1 - Plnorm), plnorm (Rlnorm, meanlog = -1, sdlog = 3, lower=F, log=T)) All.eq(log(1 - Plogis), plogis (Rlogis, location = 12, scale = 2, lower=F, log=T)) All.eq(log(1 - Pnbinom), pnbinom (Rnbinom, size = 7, prob = .01, lower=F, log=T)) All.eq(log(1 - Pnorm), pnorm (Rnorm, mean = -1, sd = 3, lower=F, log=T)) All.eq(log(1 - Ppois), ppois (Rpois, lambda = 12, lower=F, log=T)) All.eq(log(1 - Psignrank), psignrank(Rsignrank, n = 47, lower=F, log=T)) All.eq(log(1 - Pt), pt (Rt, df = 11, lower=F, log=T)) All.eq(log(1 - Pt2), pt (Rt2,df = 1.01, lower=F, log=T)) All.eq(log(1 - Punif), punif (Runif, min = .2, max = 2, lower=F, log=T)) All.eq(log(1 - Pweibull), pweibull (Rweibull, shape = 3, scale = 2, lower=F, log=T)) All.eq(log(1 - Pwilcox), pwilcox (Rwilcox, m = 13, n = 17, lower=F, log=T)) cat("Time elapsed: ", proc.time() - .ptime,"\n")