#### d|ensity
#### p|robability (cumulative)
#### q|uantile
#### r|andom number generation
####
#### Functions for ``d/p/q/r''
###-- 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))
###--- Discrete Distributions: Simple Consistency Checks pZZ = cumsum(dZZ)
all(dpois(0:5,0) == c(1, rep(0,5)))
all(dpois(0:5,0, log=TRUE) == c(0, rep(-Inf, 5)))
## Currently, just Wilcoxon [should do this for all !]
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
##--- 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*(k+1)),
pp <- sum(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(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(k, lambda=lambda, lower.tail=FALSE))
if(!(is.logical(tst3) && tst3))
cat("lambda=", format(lambda),". k =",k, " --> tst3=", tst3,"\n")
}
##--- 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 two ways:
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(tst) && tst))
cat("n=", n,"; p =",format(p),". k =",k, " --> tst=",tst,"\n")
}
}
cat("\n")
}
##--- 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-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,Inf) == 0
all(is.nan(c(pgamma(Inf,1,Inf), pgamma(Inf,Inf,1), pgamma(Inf,Inf,Inf))))
pgamma(Inf,1,xMax) == 1 && pgamma(xMax,1,Inf) == 0
##--- Normal (& Lognormal) :
## 2 Ex. from Wichura (1988); AS 241 Applied Statistics -- would be better
all.equal(qnorm(0.25),-0.6744897501960817, tol = 1e-14)
all.equal(qnorm(0.001),-3.090232306168, tol = 1e-12)
z <- rnorm(1000); all.equal(pnorm(z), 1 - pnorm(-z), tol= 1e-15)
z <- c(-Inf,Inf,NA,NaN, rt(1000, df=2))
zsml <- 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[zsml]), pnorm(z[zsml], log=TRUE))
y <- seq(-70,0, by = 10)
cbind(y, "log(pnorm(y))"= log(pnorm(y)), "pnorm(y, log=T)"= pnorm(y, log=TRUE))
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)
## The prefixes of ALL the PDQ & R functions
PDQR <- c("beta", "binom", "cauchy", "chisq", "exp", "f",
"gamma", "geom", "hyper", "lnorm", "logis", "nbinom","norm",
"pois","signrank","t","unif","weibull","wilcox")
PQonly <- c("tukey")
###===== 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) )
(Runif <- runif (n, min = .2, max = 2) )
(Rweibull <- rweibull (n, shape = 3, scale = 2) )
(Rwilcox <- rwilcox (n, m = 13, n = 17) )
(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) )
(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.equal(Rt, qt (Pt, df = 11), tol = 1e-7 )## << !!
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.equal(Rt, qt (1- Pt, df = 11, lower=F), tol = 1e-7)
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.equal(Rt, qt (log(Pt), df = 11, log=TRUE), tol = 1e-7)
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.equal(Rt, qt (log(1- Pt), df = 11, lower=F, log=T), tol = 1e-7)#ok{prec.}
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 - 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))