Nothing
tests/qgamma-ex.R
In DPQ: Density, Probability, Quantile ('DPQ') Computations
library(DPQ)
###---> Automatically find places where qgamma() is not so precise (PR#2214) :
### For PR#2214, had '1e-8' below and found quite a bit
## see /u/maechler/R/MM/NUMERICS/dpq-functions/beta-gamma-etc/qgamma-ex.R ..
## FIXME: Timing ! --- partly these matplot() partly get quite slow ~?
source(system.file(package="DPQ", "test-tools.R", mustWork=TRUE))
## => showProc.time(), ... list_() , loadList() , readRDS_() , save2RDS()
showProc.time()
(doExtras <- DPQ:::doExtras())
(sdir <- system.file("safe", package="DPQ")) ## save directory (to read from)
### Nowadays finds cases in a special region for really small p and cutoff 1e-11 :
set.seed(47)
n <- if(doExtras) 100 else 32
res <- cbind(p=1,df=1,rE=1)[-1,]
for(M in 1:(if(doExtras) 20 else 10))
for(p in runif(n)) for(df in rlnorm(n)) {
r <- 1- pchisq(qchisq(p, df),df)/p
if(abs(r) > 1e-11) res <- rbind(res, c(p,df,r))
}
### use df in U[0,1]: finds two cases with bound 1e-11
for(p in runif(n)/2) for(df in runif(n)) {
qq <- qchisq(p, df)
if(qq > 0 && p > 0) {
r <- 1- pchisq(qq, df) / p
if(abs(r) > 1e-11) res <- rbind(res, c(p,df,r))
}
}
### now df very close to 0 : ==> finds more cases
for(p in sort(c(runif(64)/2, exp(-(1+rlnorm(256))))))
for(df in 2^-rlnorm(256, mean=2, sdlog=1.5)) {
qq <- qchisq(p, df)
if(qq > 0 && p > 0) {
r <- 1- pchisq(qq, df) / p
if(abs(r) > 1e-11) res <- rbind(res, c(p,df,r))
}
}
showProc.time()
require(graphics)
if(!dev.interactive(orNone=TRUE)) pdf("qgamma-appr.pdf")
eaxis <- sfsmisc::eaxis
showProc.time()
## if(nrow(res) > 0) {
cat("Found inaccurate examples where pchisq(qchisq(p, df),df) != p\n")
## sort in p, then df:
res <- res[order(res[,"p"], res[,"df"]), ]
rE <- res[,"rE"]
if(nrow(res) > 20) { hist(rE, breaks = 30); rug(rE) }
plot(res[,1:2])##--> quite interesting : all along one curve
## p <= 1/2 and df <= 1 (about) !!
res <- cbind(res, nDig = round(-log10(abs(rE)), 1))
print(res, digits=12)
if(requireNamespace("scatterplot3d")) {
scatterplot3d::scatterplot3d(res[,1:3], type ='h') ## quite interesting:
## the inaccurate (p,df) points are on nice monotone curve !!!
## this is *less* revealing
scatterplot3d::scatterplot3d(res[,c("p","df","nDig")], type ='h')
}
rL <- res[abs(res[,'rE']) > 1e-9,]
rL <- rL[order(rL[,1],rL[,2]),]
rL
plot(rL[,1:2], type = "b", main = "inaccurate pchisq/qchisq pairs")
plot(rL[,1:2], type = "b", log = "x", ylim = range(0, rL[,"df"]),
xaxt = "n",
main = "inaccurate pchisq/qchisq pairs"); abline(h = 0, lty=2)
## aha -- a perfect line !!
lines(res[,1:2], col = adjustcolor(1, 0.5))
eaxis(1); axis(1, at = 1/2)
d <- as.data.frame(res)
plot (df ~ log(p), data = d, type = "b", cex=1/4, col="gray")
points(df ~ log(p), data = as.data.frame(rL), col=2, cex = 1/2)
summary(fm <- lm (df ~ log(p), data = d, weights = -log(abs(rE))))
## R^2 = 0.9998
p0 <- 2^seq(-50,-1, by=1/8)
dN <- data.frame(p = p0,
df = predict(fm, newdata = data.frame(p = p0)))
rE <- with(dN, 1- pchisq(qchisq(p, df),df)/p)
dN <- cbind(dN, rE = rE, nDig = round(-log10(abs(rE)), 1))
print(dN, digits=10)
## } ## only when we find inaccurate regions
showProc.time()
## Oops: another qgamma() / qchisq() problem: mostly NaN's == all solved now
curve(qgamma(x, 20), 1e-16, 1e-10, log='x')
curve(qgamma(x, 20), 1e-300, .99 , log='xy') # and add the critical region from above:
abline(v=c(1e-16, 1e-10), col="light blue")
curve(qgamma(x, 20), 1e-26, 1e-07, log='x')
##-> now using log=TRUE in same region:
curve(qgamma(x, 20, log=TRUE), -38, -16)## no problem!!
curve(qgamma(exp(x), 20), add=TRUE, col="green3", n=2001)
## had problem here, but no longer !
##--> Further fix for qgamma: when 'x' is very small: use "log=TRUE of log(x)"!
## had bug (gave NaN), but no longer:
(q_12 <- qgamma(1e-12, 20))
all.equal(1e-12, pgamma(q_12, 20), tol=0)# show rel.err (Lnx 64-bit: 4.04e-16)
stopifnot(
all.equal(1e-12, pgamma(q_12, 20), tolerance = 1e-14)
)
## --- Nice graphic : --- but amazingly *S..L..O..W*
p.qgammaSml <- function(from= 1e-110, to = 1e-5, ylim = c(0.4, 1000),
n = 201, k.lab = 3,
a1 = c(10, seq(10.1,20, by=.2), 21:105),
a2 = seq(110,330, by=10),
a3 = seq(350,1600, by=50))
{
## Purpose: nice qgamma() lines ``for small x'' aka p
## ----------------------------------------------------------------------
## Author: Martin Maechler, Date: 22 Mar 2004, 14:23
x <- exp(seq(log(from), log(to), length = n))
op <- par(las=1, lab = c(10,10, 7), xaxs = "i", mex = 0.8)
on.exit(par(op))
plot(x, qgamma(x, a1[1]), log="xy", ylim=ylim, type='l', xaxt = "n",
main = paste("qgamma(x, a) for very small x, a in [",
formatC(a1[1]),", ",formatC(max(a1,a2,a3)),"] - log-log", sep=''),
sub = R.version.string)
lab.x <- pretty(log10(c(from,to)), 20)
axis(1, at=10^lab.x, lab = paste("10^",formatC(lab.x),sep=''))
if(is.nan(qgamma(1e-12, 20)))
text(1e-60, 20, "all NaN", cex = 2)
if(!is.finite(qgamma(1e-140, 155)))
text(1e-240, 5, "all +Inf", cex = 2)
lines.txt <- function(a.s, col = par("col")) {
col <- rep(col, length=length(a.s))
for(i in seq(along=a.s)) {
qx <- qgamma(x, (a <- a.s[i]))
if(i %% k.lab == 0 &&
any(ifi <- is.finite(qx) & qx >= ylim[1])) {
ik <- (i%%(2*k.lab))/k.lab # = 0 or 1
j <- quantile(which(ifi), c(.02,(1:3)/4+ ik/10, .98))
## "segments" around the labels :
i0 <- 1
for(jj in j) {
ii <- i0:(jj-1)
i2 <- jj + -1:1
lines(x[ii], qx[ii], col=col[i])
lines(x[i2], qx[i2], col=col[i], type = 'c')
i0 <- jj+1
}
text(x[j], qx[j], formatC(a), col= "gray40", cex = 0.8)
}
else
lines(x, qx, col=col[i])
}
}
oo <- options(warn = -1)
lines.txt(a1[-1])
lines.txt(a2, col= 2)
lines.txt(a3, col= rainbow(length(a3), .8, .8,
start = (max(a3)-min(a3))/(1+max(a3))))
invisible(options(oo))
}
showProc.time()
p.qgammaSml()
p.qgammaSml(1e-300)
p.qgammaSml(1e-300,1e-50, a2= seq(100,360, by=4), a3=seq(350,1500, by=10))
showProc.time()
## The "upper" problematic corner:
p.qgammaSml(1e-19, 1e-3, a2=NULL,a3=NULL, ylim=c(.1,20))
p.qgammaSml(1e-19, 1e-3, a2=seq(1,12, by=.04), ylim=c(.1,20),a3=NULL,k.lab=10)
## now shows the problem (quite well):
## could it be in pgamma()'s inaccuracy, leading to qgamma() bias ?
aa <- c(seq(9, 22, by=0.25),seq(22.3,40,by=0.4))
caa <- formatC(range(aa))
sfsmisc::mult.fig(2)
curve(pgamma(x, sh=aa[1]), 0.5, 20, log = 'xy', ylim = c(1e-60, .2),
main = sprintf("pgamma(x, a) for a in [%s,%s]", caa[1],caa[2]))
for(sh in aa) curve(pgamma(x, sh), add = TRUE, col=2)
abline(h=c(1e-15), col="light blue", lty=2)
curve(pgamma(x, sh=aa[1]), 0.5, 20, log = 'xy', ylim = c(1e-15, .8),
main = sprintf("pgamma(x, a) for a in [%s,%s]", caa[1],caa[2]))
for(sh in aa) curve(pgamma(x, sh), add = TRUE, col=2)
## the "border curve" between "Pearson" and "Continued fraction (upper tail)"
## in pgamma.c :
curve(pgamma(max(1,x), x), add = TRUE, col=4)
## ==> pgamma() is perfect here {series expansion up to eps_C accuracy}!
aa <- c(seq(9, 22, by=0.25),seq(22.3,40.4,by=0.4))
p.qgammaSml(1e-24, 1e-5, a1=aa, a2=NULL,a3=NULL, ylim=c(.8,8))
## -------- save the above?
aa1 <- c(aa,seq(40.5,90, by=0.5))
p.qgammaSml(1e-60, 1e-5, a1=aa1, a2=NULL,a3=NULL, ylim=c(.9, 16))
aa2 <- c(aa1, seq(91,150, by= 1))
p.qgammaSml(1e-90, 1e-5, a1=aa2, a2=NULL,a3=NULL, ylim=c(.9, 35))
aa3 <- c(aa2, seq(150,250, by= 2), seq(253, 400, by=5))
p.qgammaSml(1e-200, 1e-5, a1=aa3, a2=NULL,a3=NULL, ylim=c(.9, 100))
p.qgammaSml(1e-200, 1e-5, a1=aa3, a2=NULL,a3=NULL, ylim=c(.9, 200),k.lab=9e9)
p.qgammaSml(1e-60, 1e-5, a1=aa3, a2=NULL,a3=NULL, ylim=c(.9, 200),k.lab=9e9)
showProc.time()
## lower a \> 10
curve(qgamma(x, 19), 1e-14, 1e-9, log='x')
curve(qgamma(x, 18), 1e-14, 1e-9, log='x')
curve(qgamma(x, 15), 1e-11, 5e-9, log='x')
curve(qgamma(x, 13), 5e-10, 1e-8, log='x')
curve(qgamma(x, 11), 1e-8, 5e-8, log='x')
curve(qgamma(x, 10.5), 4.2e-8, 6e-8, log='x')
curve(qgamma(x, 10.3), 6e-8, 7e-8, log='x')
curve(qgamma(x, 10.2), 7.1e-8, 7.6e-8, log='x')
curve(qgamma(x, 10.15),7.7e-8, 7.9e-8, log='x')
curve(qgamma(x, 10.14),7.88e-8,7.92e-8, log='x',n=10001)
## no more problems for smaller a!! here:
curve(qgamma(x, 10.13), 1e-10, 5e-4, log='x',n=20001)
curve(qgamma(x, 10.12), 1e-10, 5e-4, log='x',n=20001)
curve(qgamma(x, 10.1), 1e-10, 5e-4, log='x',n=20001)
showProc.time()
##--- the "+Inf" / premature "0" case:
curve(qgamma(x, 155, log=TRUE), -1500, 0, log='y', n=2001,col=2)
curve(qgamma(x, 1e3, log=TRUE), -1500, 0, log='y', n=2001,col=2)
## now works, but slowly and with kink
curve(qgamma (x, 1e5, log=TRUE), -3e5, 0, log='y', n=2001,col=2,lwd=3)
curve(qgammaAppr(x, 1e5, log=TRUE), add = TRUE, n=2001, col="blue",lwd=.4)
## --- curves are almost "identical"
## ===> the kink *does* come from the initial approx... hmm
## still "identical"
curve(qgamma (x, 1e4, log=TRUE), -3e4, 0, log='y', n=2001,col=2)
curve(qgammaAppr(x, 1e4, log=TRUE), add = TRUE, n=2001, col="tomato3")
## now see some difference (approx. has kink at ~ -165)
curve(qgamma (x, 100, log=TRUE), -200, 0, log='y', n=2001,col=2)
curve(qgammaAppr(x, 100, log=TRUE), add = TRUE, n=2001, col="tomato3")
##
(kk <- 100 * 2/1.24)# 161.29
curve(qgamma (x, 100, log=TRUE), -1.1*kk, -.95*kk, log='y', n=2001,col=2)
curve(qgammaAppr(x, 100, log=TRUE), add = TRUE, n=2001, col="tomato3")
abline(v = -kk, col='blue', lty=2)# exactly: kink is at a * 2 / 1.24 = a / .62
curve(qgammaAppr(x - 100/.62, 100,log=TRUE), -1e-3, +1e-3)
showProc.time()
p.qgammaLog <- function(alpha, xl.f = 1.5, xr.f = 0.4, n = 2001)
{
## Purpose:
## ----------------------------------------------------------------------
## Arguments:
## ----------------------------------------------------------------------
## Author: Martin Maechler, Date: 30 Mar 2004, 18:44
kk <- -alpha / .62 # = (alpha * 2) / (-1.24)
curve(qgamma(x, alpha, log=TRUE), xl.f*kk, xr.f*kk, log='y',
n=n, col=2, lwd=3.6, lty = 4,
main= paste("qgamma(x, alpha=",formatC(alpha,digits=10),", log = TRUE)"))
lines(kk, qgamma(kk, alpha, log=TRUE), type = 'h', lty = 3)
curve(qgamma (exp(x), alpha), add = TRUE, col="orange", n=n, lwd= 2)
curve(qgammaAppr(x, alpha, log=TRUE), add = TRUE, col=3, n=n,lwd = .4)
}
showProc.time()
p.qgammaLog(25)
p.qgammaLog(16)# ~ [-25, -20]
p.qgammaLog(12, 1.2, 0.8)# small problem remaining
p.qgammaLog(11, 1.2, 0.8)# even smaller
p.qgammaLog(10.5, 1.1, 0.9)# even smaller
p.qgammaLog(10.25, 1.1, 0.9)# even smaller
## 2019-08: __nothing__ visible from here on:
p.qgammaLog(10.18, 1.02, 0.98)# even smaller
p.qgammaLog(10.15, 1.02, 0.98)# even smaller
p.qgammaLog(10.14, 1.001, 0.999)# even smaller
p.qgammaLog(10.139, 1.0002, 0.9998)#
p.qgammaLog(10.138, 1.0002, 0.9998)#
p.qgammaLog(10.137, 1.00001, 0.99999)#
p.qgammaLog(10.13699, 1.0000001, 0.9999999)#
p.qgammaLog(10.1369899, 1.0000001, 0.9999999)#
p.qgammaLog(10.1369894, 1.0000001, 0.9999999)#
p.qgammaLog(10.1369893, 1.0000001, 0.9999999)# even smaller at -16.34998
showProc.time()
##-- here is the boundary --- for 64-bit AMD Opteron ---
## and for 32-bit AMD Athlon
p.qgammaLog(10.1369892, 1.0000001, 0.9999999)# no more
p.qgammaLog(10.136989, 1.0000001, 0.9999999)#
p.qgammaLog(10.136988, 1.0000001, 0.9999999)#
p.qgammaLog(10.136985, 1.0000001, 0.9999999)#
p.qgammaLog(10.13698, 1.0000001, 0.9999999)#
p.qgammaLog(10.13697, 1.0000001, 0.9999999)#
p.qgammaLog(10.13695, 1.0000001, 0.9999999)#
p.qgammaLog(10.1368, 1.000001, 0.999999)#
p.qgammaLog(10.1365, 1.000001, 0.999999)#
p.qgammaLog(10.136, 1.000001, 0.999999)#
p.qgammaLog(10.125, 1.1, 0.9)# --- see it now
p.qgammaLog(10, 1.2, 0.8)
p.qgammaLog(9)
showProc.time()
## For large alpha: show difference to see problem better
## ---> for alpha >= 10, the x problem starts *roughly* at x = -0.8*alpha
##
sfsmisc::mult.fig(2)
curve(qgammaAppr(x, 5, log=TRUE), - 8.1, -8, n=2001)
curve(qgammaAppr(x- 5/.62, 5, log=TRUE), -1e-15, 0)
## is the kink from pgamma() ? : no: this looks fine,
curve(pgamma(x, 1e5, log=TRUE), 1, 2e5, log='x', n=2001,col=2)
## and this does too:
curve( dgamma(x, 1e5), .5e5, 2e5); par(new=TRUE)
curve( dgamma(x, 1e5, log=TRUE), .5e5, 2e5, col=2, yaxt="n")
axis(4,col.axis=2); par(new=TRUE)
curve( pgamma(x, 1e5), .5e5, 2e5, n=2001, col=3); par(new=TRUE)
curve( pgamma(x, 1e5, log=TRUE), .5e5, 2e5, n=2001, col=4); par(new=TRUE)
curve(-pgamma(x, 1e5, log=TRUE,lower=FALSE), .5e5, 2e5, n=2001, col=4)
## all looking nice
x <- 10^seq(2,6, length=4001)
qx <- qgamma(pgamma(x, 1e5, log=TRUE), 1e5, log=TRUE)
plot(x, qx, type ='l', col=2, asp = 1); abline(0,1, lty=3)
showProc.time()
###------------- Approximations of qgamma() ------
##
## source("/u/maechler/R/MM/NUMERICS/dpq-functions/qchisqAppr.R")
##--> qchisqAppr()
##--> qchisqWH [ = Wilson Hilferty ]
##--> qchisqKG [ = Kennedy & Gentle's improvements "a la AS 91" ]
## dyn.load("/u/maechler/R/MM/NUMERICS/dpq-functions/qchisq_appr.so")
## Consider the two different implementations of
## lgamma1p(a) := lgamma(1+a) == log(gamma(1+a) == log(a*gamma(a)) "stable":
lseq <- if(requireNamespace("sfsmisc")) sfsmisc::lseq else
function(from, to, length) exp(seq(log(from), log(to), length.out = length))
if(require("Rmpfr")) { ##---------------- MPFR numbers -------------------------
.mpfr.all.eq <- Rmpfr::all.equal
AllEq <- function(target, current, ...)
.mpfr.all.eq(target, current, ...,
formatFUN = function(x, ...) Rmpfr::format(x, digits = 9))
print(gammaE <- Const("gamma",200)); pi. <- Const("pi",200)
print(a0 <- (gammaE^2 + pi.^2/6)/2)
print(psi2.1 <- -2*zeta(mpfr(3,200)))# == psigamma(1,2) =~ -2.4041138
print(a1 <- (psi2.1 - gammaE*(pi.^2/2 + gammaE^2))/6)
x <- lseq(1e-30, 0.8, length = if(doExtras) 1000 else 125)
x. <- mpfr(x, 200)
xct. <- log(x. * gamma(x.)) ## using MPFR arithmetic .. no overflow ...
xc2. <- log(x.) + lgamma(x.)## (ditto)
print(AllEq(xct., xc2., tol = 0)) # 3.15779......e-57
xct <- as.numeric(xct.)
stopifnot(exprs = {
AllEq(xct., xc2., tol = 1e-45)
AllEq(xct , xc2., tol = 1e-15)
##
all.equal(lgamma1p(x), lgamma1p(x, tol= 1e-16), tol=0)
## -> no difference; i.e., default tol = 1e-14 seems fine enough!
})
showProc.time()
m.appr <- cbind(log(x*gamma(x)), lgamma(1+x), log(x) + lgamma(x),
lgamma1p.(x, k=1, cut=3e-6),
lgamma1p.(x, k=2, cut=1e-4),
lgamma1p.(x, k=3, cut=8e-4),
lgamma1p(x))#, tol= 1e-14), # = default
eMat <- m.appr - xct # absolute error
## Relative errors:
str(reMat. <- m.appr/xct. - 1)
str(reMat <- as(reMat., "array")) # as(., "matrix") fails in older versions
matplot(x, eMat , log="x", type="l", lty=1) #-> problematic log(x) + lgamma(x) for "large"
matplot(x, abs( eMat), log="xy", type="l", lty=1) #-> but good for small; lgamma1p is much better
matplot(x, abs(reMat), log="xy", type="l", lty=1)
abline(v= 3.47548562941137e-08, col = "gray80", lwd=3)#<- the cutoff value of lgamma1p()
##---> should use earlier cutoff!
## zoom in:
matplot(x, abs(reMat), log="xy", type="l", col=1:7, lty=1,
lwd=2, xlim=c(8e-9, 1e-3), ylim = c(1e-18, 1e-7), axes=FALSE, frame=TRUE,
main = expression(lgamma1p(x) == log(Gamma(x+1)) ~~~ "approximations"
~~~ abs(rel.Err(.))))
eaxis(1); eaxis(2)
abline(v= 3.47548562941137e-08, col = "gray80", lwd=3)#<- the cutoff value of lgamma1p()
abline(h= c(1,2,4)*.Machine$double.eps, lty=3, col="skyblue")
legend("topright", col=1:7, lty=1,lwd=2,
c("log(x*gamma(x))", "lgamma(1+x)", "log(x) + lgamma(x)",
"lgamma1p.(x, k=1, c=3e-6)",
"lgamma1p.(x, k=2, c=1e-4)",
"lgamma1p.(x, k=3, c=8e-4)",
"lgamma1p(x)"), bty="n", ncol=2)
abline(v = c(3e-6, 1e-4, 8e-4), col=4:6, lty=2, lwd=1/2)
## FIXME: do the same for the lgaamma1p_series()
## rm(x., xct., xc2., reMat., eMat, AllEq)
detach("package:Rmpfr")
showProc.time()
} ## if( MPFR ) ----------------------------------------------------------------
## ../R/qchisqAppr.R -- talks about the "small shape" qgamma() approximation
## ----------------- --> .qgammaApprBnd() :
curve(.qgammaApprBnd, 1e-18, 1e-15, col=2)
abline(h=0, col="gray70", lty=2)
eps.c <- .Machine$double.eps
axis(3,at=(1:3)* eps.c,
label=expression(epsilon[c], 2*epsilon[c], 3*epsilon[c]))
(rt.b <- uniroot(.qgammaApprBnd, c(1,3)*eps.c, tol=1e-12))
rt.b$root ## 3.954775e-16
rt.b$root / eps.c ## 1.781072
##==> for a < 1.781*eps, bnd > 0 ==> we have log(p) < bnd for all p
## otherwise, we should effectively 'test'
curve(.qgammaApprBnd, 1e-16, 1e-10, log="x", col=2)
showProc.time()
## source ("/u/maechler/R/MM/NUMERICS/dpq-functions/beta-gamma-etc/qgamma-fn.R")
## ##--> qchisqAppr.R() -- which has 'kind = ' argument!
## ##--> qgamma.R()
p.qchi.appr <-
function(x, qm= { m <- cbind(qchisq(x, df, log=TRUE),
sapply(knds, function(kind)
qchisqAppr.R(x,df,log=TRUE,kind=kind)))
colnames(m) <- c("True", "default", knds[-1])
m },
df,
knds = list(NULL,"chi.small", "WH", "p1WH", "df.small"),
call = match.call(), main = deparse(call), log = "y", ...)
{
## Purpose:
## ----------------------------------------------------------------------
## Arguments:
## ----------------------------------------------------------------------
## Author: Martin Maechler, Date: 25 Mar 2004, 22:08
col <- c(2,1,3:6)
lty <- c(1,3,1,1,1,1)
lwd <- c(2,2,1,1,1,1)
matplot(x, qm, col=col, lty=lty, lwd=lwd, log = log,
type = 'l', main = main, ...)
y0 <- c( .02, .98) %*% par('usr')[3:4]
if(par("ylog")) y0 <- 10^y0
legend(min(x), y0, colnames(qm), col=col, lty=lty, lwd=lwd)
invisible(list(x=x, qm=qm, call=match.call()))
}
pD.chi.appr <- function(pqr, err.kind=c("relative", "absolute"),
type = "l", log = "y",
lwds = c(2, rep(1, k-1)),
cols = seq(along=lwds),
ltys = rep(1,k),
...)
{
## Purpose: Plot Difference from "True" qchisq()
## ----------------------------------------------------------------------
## Arguments: pqr: a list as resulting from p.chi.appr()
## ----------------------------------------------------------------------
## Author: Martin Maechler, Date: 31 Mar 2004, 09:38
err.kind <- match.arg(err.kind)
if(!is.list(pqr) || !is.numeric(k <- ncol(pqr$qm)-1) || k <= 0)
stop("invalid first argument 'pqr'")
with(pqr, {
err <- abs(if(err.kind == "relative")
(1- qm[,-1] / qm[,1]) else (qm[,-1] - qm[,1]))
matplot(x, err, ylab = "",
main = paste(err.kind,"Error from", deparse(call)),
type= type, log= log, lty=ltys, col=cols, lwd=lwds, ...)
legend(par("xaxp")[1], par("yaxp")[2], colnames(qm)[-1],
lty=ltys, col=cols, lwd=lwds)
})
}
## if(FALSE) # you can manually set
## do.pqchi <- TRUE # before source()ing this file
## if(!exists("do.pqchi") || !is.logical(do.pqchi))
## do.pqchi <- interactive()
## if(do.pqchi) { #------------- FIXME look at speed up or cache indeed !?? <<<<<
## pqFile <- "/u/maechler/R/MM/NUMERICS/dpq-functions/pqchi.rda"
## ## ls -l ................... 1325446 Nov 2 2009 pqchi.rda
## if(file.exists(pqFile)) {
## attach(pqFile) ## it loads more than we create here __FIXME__
## print(ls(2, all.names=TRUE))
## ## [1] "pq.1" "pq.25" "pq.25." "pq.31" "pq.33" "pq.33." "pq.33.2"
## ## [8] "pq.33.3" "pq.33.4" "pq.5" "pq.5." "pq1" "pq1." "pq1.95"
## ## [15] "pq1.95." "pq1.95.2" "pq10" "pq10." "pq10.2" "pq100" "pq2"
## ## [22] "pq2." "pq2.05" "pq2.05." "pq2.05.2" "pq2.5" "pq2.5." "pq2.5.2"
## ## [29] "pq200" "pq2L" "pq4" "pq4." "pq4.2" "pqFile"
## }
showProc.time()
x <- seq(-300, -.01, length=501)
all.equal(qchisqAppr.R(x, 200, log=TRUE),
qchisqAppr (x, 200, log=TRUE),tol=0)
## 4.48 e-16 / TRUE (Opteron)
all.equal(qchisqAppr.R(x, 2, log=TRUE),
qchisqAppr (x, 2, log=TRUE),tol=0)
## 3.90 e-16 / TRUE (Opteron)
all.equal(qchisqAppr.R(x, 0.1, log=TRUE),
qchisqAppr (x, 0.1, log=TRUE),tol=0)
## 7.15 e-15 / 1.179e-8 !!!!! (Opteron)
pq200 <- p.qchi.appr(x = seq(-300, -.01, length=501), df = 200)
pq100 <- p.qchi.appr(x = seq(-160, -.01, length=501), df = 100)
## after (slow) computing, quickly repeat:
with(pq200, p.qchi.appr(x=x, qm=qm, call=call))
with(pq100, p.qchi.appr(x=x, qm=qm, call=call))
## this "hangs forever" -- before I introduced 'maxit' (for 'nu.small'):
pq10 <- p.qchi.appr(x = seq(-12, -.005, length=501), df = 10)
## want to see the jump:
pq10. <- p.qchi.appr(x = seq(-10, -4, length=501), df = 10)
pq10.2<- p.qchi.appr(x = seq(-8.5,-7.5, length=501), df = 10)
with(pq10.2, p.qchi.appr(x=x, qm=qm, call=call))
pq2.5 <- p.qchi.appr(x = seq(-3.4, -.01, length=501), df = 2.5)
pq2.5. <- p.qchi.appr(x = seq(-2.1, -1.8, length=901), df = 2.5)#the jump
## what about p1WH (which is fantastic for df=2)?
pq2.5.2<- p.qchi.appr(x = seq(-0.5, -1e-3, length=901), df = 2.5)
with(pq2.5, p.qchi.appr(x=x, qm=qm, call=call))
with(pq2.5.2, p.qchi.appr(x=x, qm=qm, call=call))
pD.chi.appr(pq2.5.2)# nothing special
pq2.05 <- p.qchi.appr(x = seq(-3.4, -.01, length=501), df = 2.05)
pq2.05. <- p.qchi.appr(x = seq(-2.5, -1.5, length=901), df = 2.05)#the jump
## ^^ the jump from chi.small to WH is much too late here
## what about p1WH (which is fantastic for df=2)?
pq2.05.2<- p.qchi.appr(x = seq(-0.4, -1e-5, length=901), df = 2.05)
pD.chi.appr(pq2.05.2) # p1WH is starting to become better
# and the jump (WH -> p1WH) is too late
with(pq2.05, p.qchi.appr(x=x, qm=qm, call=call))
with(pq2.05.2, p.qchi.appr(x=x, qm=qm, call=call))
## Here, all are 'ok' (but "nu.small"):
pq2L <- p.qchi.appr(seq(-300, -.01, length=201), df = 2)
pq2 <- p.qchi.appr(x = seq(-5, -.001, length=501), df = 2)
pq2. <- p.qchi.appr(x = seq(-2.5, -1, length=901), df = 2)
with(pq2., p.qchi.appr(x=x, qm=qm, call=call))
pq4 <- p.qchi.appr(x = seq(-8, -0.01, length = 501), df = 4)
summary(warnings())
pq4.2 <- p.qchi.appr(x = seq(-0.1, -1e-04, length = 901), df = 4)
pq1.95 <- p.qchi.appr(x = seq(-3., -.01, length=501), df = 1.95)
pq1.95. <- p.qchi.appr(x = seq(-2.2, -1.5, length=901), df = 1.95)#the jump -1.57
## ^^ the jump from chi.small to WH is *much* too late here
## what about p1WH (which is fantastic for df=2)?
pq1.95.2<- p.qchi.appr(x = seq(-0.4, -1e-7, length=901), df = 1.95)
pD.chi.appr(pq1.95.2) # p1WH is starting to become better
# and the jump (WH -> p1WH) is too late
with(pq1.95, p.qchi.appr(x=x, qm=qm, call=call))
with(pq1.95.2, p.qchi.appr(x=x, qm=qm, call=call))
pq1 <- p.qchi.appr(x = seq(-4, -.001, length=501), df = 1)
pq1. <- p.qchi.appr(x = seq(-1.8, -.5, length=901), df = 1)
with(pq1., p.qchi.appr(x=x, qm=qm, call=call))
pq.5 <- p.qchi.appr(x = seq(-1.5, -.001, length=501), df = 0.5)
pq.5. <- p.qchi.appr(x = seq(-0.8, -0.2, length=901), df = 0.5, ylim=c(.04,.6))
pq.33 <- p.qchi.appr(x= seq(-0.9, -.001,length=501), df= 0.33)
pq.33. <- p.qchi.appr(x= seq(-0.4, -0.02,length=901), df= 0.33)
pq.33.2<- p.qchi.appr(x= seq(-0.4, -0.2, length=901), df= 0.33, ylim=c(.15,.60))
with(pq.33.2, p.qchi.appr(x=x, qm=qm, call=call,ylim=c(.15,.45)))
pq.33.3<- p.qchi.appr(x= seq(-0.4, -0.005, length=901), df= 0.33, ylim=c(.15, 4.00))
with(pq.33.3, p.qchi.appr(x=x, qm=qm, call=call))#,ylim=c(.15, 8)))
pq.33.4<- p.qchi.appr(x= seq(-0.003, -1e-6, length=901), df= 0.33,ylim=c(5,25))
with(pq.33.4, p.qchi.appr(x=x, qm=qm, call=call,ylim=c(5,25)))
## nu <= 0.32 is the "magic" border of Best & Roberts
pq.31 <- p.qchi.appr(x = seq(-0.45,-.010, length=501), df = 0.31)
with(pq.31, p.qchi.appr(x=x, qm=qm, call=call))
pq.25 <- p.qchi.appr(x = seq(-0.3, -0.02, length=901), df = 0.25)
pq.1 <- p.qchi.appr(x = seq(-0.16, -0.01, length=901), df = 0.1)
with(pq.1, p.qchi.appr(x=x, qm=qm, call=call))
showProc.time()
## if(!file.exists(pqFile)) # don't overwrite for now (as it contains pq2L ,
## save(list=ls(pat="^pq"), file = pqFile)
##}## end if(do.pqchi){ only if interactive } ======================================
pD.chi.appr(pq2L, "abs")
pD.chi.appr(pq2L, "rel")
## --> want only much smaller x-range:
pD.chi.appr(pq2,"abs")#--> fantastic p1WH
pD.chi.appr(pq2) # (ditto)
pD.chi.appr(pq4.2)# p1WH: only at very right
pD.chi.appr(pq4.2, xlim=c(-.016,0))# p1WH: only at very right
## no Newton step here, eg:
(qgA100 <- qgammaAppr(1e-100, 100))
(qg.100 <- qgamma (1e-100, 100))
all.equal(qgA100, qg.100)
## too much different
dgamma(1e-100, 100, log = TRUE)# -23154.7 i.e., "non-log" is 0
qgamma.R(1e-100, 100, verbose = TRUE)#-> final Newton fails!
## but here, the final Newton iterations do work :
x <- qgamma.R(log(1e-100), 100, log = TRUE, verbose = TRUE)
pgamma(x, 100) # = 1e-100 ! perfect
showProc.time()
###--> Use this to devise an improved final Newton algorithm !!!
## From: Prof Brian Ripley <ripley@stats.ox.ac.uk>
## To: skylab.gupta@gmail.com
## Cc: R-bugs@biostat.ku.dk, r-devel@stat.math.ethz.ch
## Subject: Re: [Rd] qgamma inaccuracy (PR#12324)
## Date: Tue, 12 Aug 2008 20:50:50 +0100 (BST)
## This is a really extreme usage. AFAICS the code works well enough down to
## shape=1e-10 or so, e.g.
qgamma(1e-10, 5e-11, lower.tail=FALSE)
## [1] 0.08237203
## in R 2.9.. .. 2.10.0 -- with an accuracy warning {which is *wrong*!}
## I would be interested to know what substantive problem you were trying to
## solve here that required such values.
## I am pretty sure that a completely different algorithm will be required.
## MM: It looks like this is basis for a new algo:
a <- 1e-14
gammaE <- 0.57721566490153286060651209008240243079 # Euler's gamma
curve(pgamma(x, a, lower.tail=FALSE)/a + log(x) + gammaE, 1e-300, 1e-1, log="",n=1000)
curve(pgamma(x, a, lower.tail=FALSE)/a + log(x) + gammaE - x, 1e-300, 1e-1, log="",n=1000)
## ==> Q = 1 - P = pgamma(x,a, lower=FALSE) ~= a*(-log(x) - gammaE + x - x^2/4)
## i.e., Q ~= -a(log(x) + gammaE { -x + x^2/4 }
## -Q/a - gammaE ~= log(x) { -x + x^2/4 }
## ==> x ~= exp(- (Q/a + gammaE))
## e.g., example below:
Q <- 1e-100; a <- 5e-101
## MM: Find inverse :
str(r.a <- uniroot(function(x) pgamma(x,a,lower.tail=FALSE) - Q,
int = c(0.01, 0.1), tol=1e-20))
dput(x0 <- r.a$root) ## 0.0823720296207203
(x1 <- exp(- (Q/a + gammaE)))## 0.07598528 .. not so good
qgammaApprSmallP(Q, a, lower.tail=FALSE)## ~= 0.07598528 -- the same!!
pgamma(x0, a, lower.tail=FALSE) ## 1.00000e-100.
pgamma(x1, a, lower.tail=FALSE) ## 1.03728e-100 ... hmm "close"
##
## MM: -- now look at the bigger picture
p.qg.2a <- function(l2x.min= -15, l2x.max = -100, n = 501,
do.offset = FALSE,
type = "o", log = "x", cex = 0.6, ...) {
x.log <- any("x" == strsplit(log,"")[[1]])
x <- if(x.log) 2^seq(l2x.min, l2x.max, length=n)
else seq(2^l2x.min,2^l2x.max, length=n)
if(do.offset)
plot(x, qgamma(2*x, x, lower.tail=FALSE) - 0.0823720296206873,
type=type, cex=cex, log=log, ...)
else plot(x, qgamma(2*x, x, lower.tail=FALSE),
type=type, cex=cex, log=log, ...)
}
p.qg.2a() # was "very bad" in R <= 2.10.0, now --> 0.082372... "perfect" smooth
## still a little ---but acceptable--- remaining inaccuracy ...zooming in:
p.qg.2a(-43,-55, do.offset=TRUE)
p.qg.2a(,-1024)
p.qg.2a(,-1024, log="", pch=".")## linear in x !!
## zoom in at the limit
p.qg.2a(-30,-1024, do.offset=TRUE, ylim = 1e-11*c(-1,1))
p.qg.2a(-33,-1024, do.offset=TRUE, ylim = 1e-12*c(-1,1))
p.qg.2a(-33,-1024, do.offset=TRUE, ylim = 1e-13*c(-1,1))
a <- 2^-(7:900)
qg <- qgamma(2*a, a, lower.tail=FALSE)
re <- 1-pgamma(qg, a, lower.tail=FALSE)/(2*a)
plot(a, re, log="x", type="b", col=2)
stopifnot(abs(re) < 2e-12) # but really, *should be a bit better
showProc.time()
## For completeness we may write that in due course, but for now (R 2.7.2) I
## suggest just issuing a warning for miniscule 'shape'.
## On Thu, 7 Aug 2008, skylab.gupta@gmail.com wrote:
## > Full_Name:
## > Version: 2.7.1 (2008-06-23)
## > OS: windows vista
## > Submission from: (NULL) (216.82.144.137)
## >
## >
## > Hello,
## >
## > I have been working with various probability distributions in R, and it seems
## > the gamma distribution is inaccurate for some inputs.
## > For example, qgamma(1e-100, 5e-101, lower.tail=FALSE) gives: 1.0. However, it
## > seems this is incorrect; I think the correct answer should be
## > 0.082372029620717283. When I check these numbers using pgamma, I get:
(qg <- qgamma(1e-100, 5e-101, lower.tail=FALSE))# 1 (wrong, originally)
## 0.08237203 now (2009-11-04), i.e. ok
pgamma(qg, 5e-101, lower.tail=FALSE)# now -> 1e-100 : ok
pgamma(0.082372029620717283,5e-101, lower.tail=FALSE)
## 1.0000000000000166e-100.
RE.pqgamma <- function(p, shape, lower.tail = TRUE, log.p = FALSE) {
## Relative Error of pgamma(qgamma(*), ..):
1 - pgamma(qgamma(p, shape, lower.tail=lower.tail, log.p=log.p),
shape=shape, lower.tail=lower.tail, log.p=log.p) / p
}
RE.qpgamma <- function(q, shape, lower.tail = TRUE, log.p = FALSE) {
## Relative Error of qgamma(pgamma(*), ..):
1 - qgamma(pgamma(q, shape, lower.tail=lower.tail, log.p=log.p),
shape=shape, lower.tail=lower.tail, log.p=log.p) / q
}
## Ok, how extreme can we get -- let a := alpha := shape --> 0 :
x <- 1e-100
a <- 2^-(7:300)# is still "ok":
plot(a, (re <- RE.pqgamma(x, a, lower.tail=FALSE)), type="b", col=2, log="x")# oops!
a <- 2^-(7:400)
plot(a, (re <- RE.pqgamma(x, a, lower.tail=FALSE)), type="b", col=2, log="x")# oops!
## Oops!
## but, it is clear
qgamma(x, 2^-400, lower.tail=FALSE)## is exactly 0
## -> it goes to 0 quickly . .. zooming in:
curve(qgamma(1e-100, x, lower.tail=FALSE), 1e-110, 1e-70, log="xy", col=2, axes=FALSE)
eaxis(1);eaxis(2)
## from when on is it exactly 0:
uniroot(function(u) qgamma(1e-100, 2^u,lower.tail=FALSE)-1e-315, c(-400, -300))$root
## -341.6941
## => use
a <- 2^-(7:341)
plot(a, (re <- RE.pqgamma(x, a, lower.tail=FALSE)),
type="l", col=2, log="x")# small glips
## zoom in:
curve(abs(RE.pqgamma(1e-100, x, lower.tail=FALSE)), 2^-341, 1e-90, log="xy", n=2000)
curve(RE.pqgamma(1e-100, x, lower.tail=FALSE), 1e-100, 10e-100, n=2000)
curve(RE.pqgamma(1e-100, x, lower.tail=FALSE), 4e-100, 6e-100, n=2000)
## Ok: at least here is a problem
RE.pqgamma(1e-100, 5e-100, lower.tail=FALSE)# -0.1538
## more general
curve(RE.pqgamma(x, 5*x, lower.tail=FALSE), 1e-100, 1e-20, n=10000, log="x")
## problem *everywhere* , starting quite early: (lesser problem at ~ 1e-16 !)
curve(RE.pqgamma(x, 5*x, lower.tail=FALSE), 1e-25, 1e-15, n=1000, log="x")
curve(RE.pqgamma(x, 5*x, lower.tail=FALSE), 1e-21, 10e-21,n=1000)
curve(RE.pqgamma(x, 5*x, lower.tail=FALSE), 2e-21, 6e-21, n=1000)
curve(RE.pqgamma(x, 5*x, lower.tail=FALSE), 4e-21, 4.5e-21, n=1000)
## and indeed, it's qgamma() that jumps here
curve(qgamma(x, 5*x, lower.tail=FALSE), 4e-21, 4.5e-21, ylim=c(.97, 1.1))
## well, looking at pgamma(), finally reveals the buglet is there first:
## There's a jump at x = 1 !!
curve(pgamma(x, 1e-30, lower=FALSE), .9999, 1.0001)
curve(pgamma(x, 1e-20, lower=FALSE), .9999, 1.0001)
curve(pgamma(x, 1e-17, lower=FALSE), .9999, 1.0001)
curve(pgamma(x, 1e-15, lower=FALSE), .9999, 1.0001)
curve(pgamma(x, 1e-13, lower=FALSE), .9999, 1.0001)
curve(pgamma(x, 1e-12, lower=FALSE), .9999, 1.0001)# still clearly visible
curve(pgamma(x, 1e-11, lower=FALSE), .9999, 1.0001)# barely visible
## for larger alpha == shape, must zoom in more and more:
curve(pgamma(x, 1e-11, lower=FALSE), .99999, 1.00001)#
curve(pgamma(x, 1e-10, lower=FALSE), .999999, 1.000001)
curve(pgamma(x, 1e-9, lower=FALSE), .9999999, 1.0000001)
curve(pgamma(x, 1e-8, lower=FALSE), .99999999, 1.00000001)
curve(pgamma(x, 1e-7, lower=FALSE), .999999999, 1.000000001)
curve(pgamma(x, 1e-6, lower=FALSE), .9999999999, 1.0000000001)
curve(pgamma(x, 1e-5, lower=FALSE), .99999999999, 1.00000000001)
curve(pgamma(x, 1e-4, lower=FALSE), .999999999999, 1.000000000001)
## now we get close to noise level:
curve(pgamma(x, 1e-3, lower=FALSE), .9999999999999, 1.0000000000001)
curve(pgamma(x, 1e-2, lower=FALSE), .99999999999999, 1.00000000000001)
showProc.time()
del.pgamma <- function(a, eps = 1e-13)
{
## Purpose: *relative* jump size at x = 1 of pgamma(x, a, lower=FALSE)
## ----------------------------------------------------------------------
## Author: Martin Maechler, Date: 5 Nov 2009, 16:08
stopifnot(a > 0, length(a) == 1, eps > 0, length(eps) == 1)
pp <- pgamma(1+c(-eps,eps), a, lower.tail = FALSE)
(pp[2] - pp[1])*2/(pp[2] + pp[1])
}
a <- lseq(1e-300, 1e-3, length=400)
dpa <- sapply(a, del.pgamma)
plot(a, -dpa, log="xy", type="l", col=2, yaxt="n");eaxis(2)
## ok, it remains constant all the way to 1e-300
## --> focus
a <- lseq(1e-40, 1e-5, length=400)
dpa <- sapply(a, del.pgamma)
plot(a, -dpa, log="xy", type="l", col=2, axes=FALSE)
eaxis(1, at = 10^-seq(5,40, by=5))
eaxis(2)
xm <- .Machine$double.xmin
pgamma(xm, shape= 1e-20)# is "practically 1" --> *most* qgamma() will be exactly 0
## how close to 1 ? ---> use upper_tail [possibly on log scale:]
pgamma(xm, shape= 1e-20, lower.tail=FALSE) # 7.078192e-18
pgamma(xm, shape= 1e-20, lower.tail=FALSE, log=TRUE) # -39.48951
## Where is the 'boundary' (from which on qgamma() must return 0, since it can't give
## xm = 2.2e-308 ) ?
a <- 2^-(7:1030)
plot(a, (p <- pgamma(xm, a, lower.tail=FALSE, log=TRUE)),
cex=.5,type ="l", col=2, log="x")
summary(lm. <- lm(p ~ log(a) + a + I(a^2))) ## coeff. of log(a) *is* 1
summary(lm. <- lm(p ~ offset(log(a)) + a + I(a^2)))
p.l <- pgamma(xm, a, lower.tail=FALSE, log=TRUE) - log(a)
dput(mean(tail(p.l))) ## 6.56218869790132
##=> pgamma(xm, a) ~= log(a) + 6.5621886979 +
## ok, to get this better, now need different a:
al <- seq(1e-300, 1e-3, length=200)
plot(al, (pl <- pgamma(xm, al, lower.tail=FALSE, log=TRUE)) - (log(al)+6.56218869790132),
cex=.5,type ="l", col=2)
summary(lm.. <- lm(pl ~ offset(log(al) + 6.56218869790132) + 0 + al + I(al^2)))
confint(lm..)
coef(lm..)
## al I(al^2)
## -353.8587 20750.8205
##=> pgamma(xm, a) ~= log(a) + 6.5621886979 - 353.858745 * a
## > Similarly, for example: -- now (2009-11-04) ok
qgamma(1e-100,0.005,lower.tail=FALSE) # = 109.36757177 instead of 219.5937..
pgamma(109.36757177007101, 0.005, lower.tail=FALSE)# = 1.4787306506694758e-52.
## > This looks completely wrong. The correct value, I think, should be
## > 219.59373661415756. In fact,
pgamma(219.59373661415756, 0.005, lower.tail=FALSE)# = 9.9999999999999558e-101.
## >
## > In fact, when I do the following in R, the results are completely wrong,
## >
a <- 5*10^-(1:40)
z1 <- qgamma(1e-100,a,lower.tail=FALSE)
(y <- pgamma(z1,a,lower.tail=FALSE))
## The value of y that I get should be close to 1e-100, but they are not:
## [1] 1.000000e-100 1.871683e-51 1.478731e-52 1.444034e-53 1.440606e-54
## [6] 1.440264e-55 1.440230e-56 1.440226e-57 1.440226e-58 1.440226e-59
summary(abs(1 - y/1e-100))
plot(a, abs(1 - y/1e-100), log="xy", type="b")
stopifnot(abs(1 - y/1e-100) < 2e-13)# max (32b Linux, P.M) = 4.186e-14
## > The correct values of z1 should be:
z1true <- c(226.97154111939946, 222.15218724493326, 219.59373661415756,
217.27485383840451, 214.98015408183574, 212.68797118872064,
210.39614286838227, 208.10445550564617, 205.81289009100664,
203.52144711679352)
all.equal(z1[1:10], z1true, tol=0)# 1.307e-16 on 32-bit (Pentium M); 1.686e-16 (64b Lnx, 2019)
stopifnot(all.equal(z1[1:10], z1true, tol = 1e-15))
showProc.time()
##>
##> With these values of z1true, we get the expected values:
(y <- pgamma(z1true, a, lower.tail=FALSE))
## [1] 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100
## > I am using the precompiled binary version of R, under Windows Vista.
## > -----------
## >> version
## > _
## > platform i386-pc-mingw32
## > arch i386
## > os mingw32
## > system i386, mingw32
## > status
## > major 2
## > minor 7.1
## > year 2008
## > month 06
## > day 23
## > svn rev 45970
## > language R
## > version.string R version 2.7.1 (2008-06-23)
## > ------------
## >
## > So, it seems qgamma is inaccurate for small probability values in the upper
## > tail, when the shape parameter is also small.
###_-- MM: Still wrong:
(xm <- 2^-1074.9999) # is less than .Machine $ double.xmin == the really x > 0
pgamma(xm, .00001)# 0.992589
qgamma(.99, .00001)##--> NaN -- should give 0 or "xmin" or so
## FIXME -- ok, now
## but
curve(qgamma(x, .001, lower=FALSE), .4, .8, n=1001, log="y")
curve(qgamma(x, 1e-5, lower=FALSE), .002, .2, n=1001, log="xy")
curve(qgamma(x, 1e-7, lower=FALSE), 1e-5, .04, n=1001, log="xy")
curve(qgamma(x, 1e-12, lower=FALSE), 1e-12, 1e-2, n=1001, log="xy")
## or
curve(qgamma(x, 1e-121, lower=FALSE), 7e-119, 8e-119,
n=2001, log="y", yaxt="n")
try( # reveals eaxis() bug ? -- for the *subnormal* numbers
eaxis(2, at = 10^-seq(304,324, by=2))
)
curve(qgamma(x, .001, lower=FALSE), .4, .6, n=1001, log="y")
curve(qgamma(x, .001, lower=FALSE), .5, .55, n=1001, log="y")
## gives a *warning* from axis() because of subnormal y-range {was error, fixed in R-devel (2018-08)}
curve(qgamma(x, .001, lower=FALSE), .52, .53, n=1001, log="y")
curve(qgamma(x, .001, lower=FALSE)*1e100, .522, .526, n=1001, log="y")
showProc.time()
Try the DPQ package in your browser
Any scripts or data that you put into this service are public.
DPQ documentation built on Dec. 5, 2023, 3:05 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
library(DPQ)
###---> Automatically find places where qgamma() is not so precise (PR#2214) :
### For PR#2214, had '1e-8' below and found quite a bit
## see /u/maechler/R/MM/NUMERICS/dpq-functions/beta-gamma-etc/qgamma-ex.R ..
## FIXME: Timing ! --- partly these matplot() partly get quite slow ~?
source(system.file(package="DPQ", "test-tools.R", mustWork=TRUE))
## => showProc.time(), ... list_() , loadList() , readRDS_() , save2RDS()
showProc.time()
(doExtras <- DPQ:::doExtras())
(sdir <- system.file("safe", package="DPQ")) ## save directory (to read from)
### Nowadays finds cases in a special region for really small p and cutoff 1e-11 :
set.seed(47)
n <- if(doExtras) 100 else 32
res <- cbind(p=1,df=1,rE=1)[-1,]
for(M in 1:(if(doExtras) 20 else 10))
for(p in runif(n)) for(df in rlnorm(n)) {
r <- 1- pchisq(qchisq(p, df),df)/p
if(abs(r) > 1e-11) res <- rbind(res, c(p,df,r))
}
### use df in U[0,1]: finds two cases with bound 1e-11
for(p in runif(n)/2) for(df in runif(n)) {
qq <- qchisq(p, df)
if(qq > 0 && p > 0) {
r <- 1- pchisq(qq, df) / p
if(abs(r) > 1e-11) res <- rbind(res, c(p,df,r))
}
}
### now df very close to 0 : ==> finds more cases
for(p in sort(c(runif(64)/2, exp(-(1+rlnorm(256))))))
for(df in 2^-rlnorm(256, mean=2, sdlog=1.5)) {
qq <- qchisq(p, df)
if(qq > 0 && p > 0) {
r <- 1- pchisq(qq, df) / p
if(abs(r) > 1e-11) res <- rbind(res, c(p,df,r))
}
}
showProc.time()
require(graphics)
if(!dev.interactive(orNone=TRUE)) pdf("qgamma-appr.pdf")
eaxis <- sfsmisc::eaxis
showProc.time()
## if(nrow(res) > 0) {
cat("Found inaccurate examples where pchisq(qchisq(p, df),df) != p\n")
## sort in p, then df:
res <- res[order(res[,"p"], res[,"df"]), ]
rE <- res[,"rE"]
if(nrow(res) > 20) { hist(rE, breaks = 30); rug(rE) }
plot(res[,1:2])##--> quite interesting : all along one curve
## p <= 1/2 and df <= 1 (about) !!
res <- cbind(res, nDig = round(-log10(abs(rE)), 1))
print(res, digits=12)
if(requireNamespace("scatterplot3d")) {
scatterplot3d::scatterplot3d(res[,1:3], type ='h') ## quite interesting:
## the inaccurate (p,df) points are on nice monotone curve !!!
## this is *less* revealing
scatterplot3d::scatterplot3d(res[,c("p","df","nDig")], type ='h')
}
rL <- res[abs(res[,'rE']) > 1e-9,]
rL <- rL[order(rL[,1],rL[,2]),]
rL
plot(rL[,1:2], type = "b", main = "inaccurate pchisq/qchisq pairs")
plot(rL[,1:2], type = "b", log = "x", ylim = range(0, rL[,"df"]),
xaxt = "n",
main = "inaccurate pchisq/qchisq pairs"); abline(h = 0, lty=2)
## aha -- a perfect line !!
lines(res[,1:2], col = adjustcolor(1, 0.5))
eaxis(1); axis(1, at = 1/2)
d <- as.data.frame(res)
plot (df ~ log(p), data = d, type = "b", cex=1/4, col="gray")
points(df ~ log(p), data = as.data.frame(rL), col=2, cex = 1/2)
summary(fm <- lm (df ~ log(p), data = d, weights = -log(abs(rE))))
## R^2 = 0.9998
p0 <- 2^seq(-50,-1, by=1/8)
dN <- data.frame(p = p0,
df = predict(fm, newdata = data.frame(p = p0)))
rE <- with(dN, 1- pchisq(qchisq(p, df),df)/p)
dN <- cbind(dN, rE = rE, nDig = round(-log10(abs(rE)), 1))
print(dN, digits=10)
## } ## only when we find inaccurate regions
showProc.time()
## Oops: another qgamma() / qchisq() problem: mostly NaN's == all solved now
curve(qgamma(x, 20), 1e-16, 1e-10, log='x')
curve(qgamma(x, 20), 1e-300, .99 , log='xy') # and add the critical region from above:
abline(v=c(1e-16, 1e-10), col="light blue")
curve(qgamma(x, 20), 1e-26, 1e-07, log='x')
##-> now using log=TRUE in same region:
curve(qgamma(x, 20, log=TRUE), -38, -16)## no problem!!
curve(qgamma(exp(x), 20), add=TRUE, col="green3", n=2001)
## had problem here, but no longer !
##--> Further fix for qgamma: when 'x' is very small: use "log=TRUE of log(x)"!
## had bug (gave NaN), but no longer:
(q_12 <- qgamma(1e-12, 20))
all.equal(1e-12, pgamma(q_12, 20), tol=0)# show rel.err (Lnx 64-bit: 4.04e-16)
stopifnot(
all.equal(1e-12, pgamma(q_12, 20), tolerance = 1e-14)
)
## --- Nice graphic : --- but amazingly *S..L..O..W*
p.qgammaSml <- function(from= 1e-110, to = 1e-5, ylim = c(0.4, 1000),
n = 201, k.lab = 3,
a1 = c(10, seq(10.1,20, by=.2), 21:105),
a2 = seq(110,330, by=10),
a3 = seq(350,1600, by=50))
{
## Purpose: nice qgamma() lines ``for small x'' aka p
## ----------------------------------------------------------------------
## Author: Martin Maechler, Date: 22 Mar 2004, 14:23
x <- exp(seq(log(from), log(to), length = n))
op <- par(las=1, lab = c(10,10, 7), xaxs = "i", mex = 0.8)
on.exit(par(op))
plot(x, qgamma(x, a1[1]), log="xy", ylim=ylim, type='l', xaxt = "n",
main = paste("qgamma(x, a) for very small x, a in [",
formatC(a1[1]),", ",formatC(max(a1,a2,a3)),"] - log-log", sep=''),
sub = R.version.string)
lab.x <- pretty(log10(c(from,to)), 20)
axis(1, at=10^lab.x, lab = paste("10^",formatC(lab.x),sep=''))
if(is.nan(qgamma(1e-12, 20)))
text(1e-60, 20, "all NaN", cex = 2)
if(!is.finite(qgamma(1e-140, 155)))
text(1e-240, 5, "all +Inf", cex = 2)
lines.txt <- function(a.s, col = par("col")) {
col <- rep(col, length=length(a.s))
for(i in seq(along=a.s)) {
qx <- qgamma(x, (a <- a.s[i]))
if(i %% k.lab == 0 &&
any(ifi <- is.finite(qx) & qx >= ylim[1])) {
ik <- (i%%(2*k.lab))/k.lab # = 0 or 1
j <- quantile(which(ifi), c(.02,(1:3)/4+ ik/10, .98))
## "segments" around the labels :
i0 <- 1
for(jj in j) {
ii <- i0:(jj-1)
i2 <- jj + -1:1
lines(x[ii], qx[ii], col=col[i])
lines(x[i2], qx[i2], col=col[i], type = 'c')
i0 <- jj+1
}
text(x[j], qx[j], formatC(a), col= "gray40", cex = 0.8)
}
else
lines(x, qx, col=col[i])
}
}
oo <- options(warn = -1)
lines.txt(a1[-1])
lines.txt(a2, col= 2)
lines.txt(a3, col= rainbow(length(a3), .8, .8,
start = (max(a3)-min(a3))/(1+max(a3))))
invisible(options(oo))
}
showProc.time()
p.qgammaSml()
p.qgammaSml(1e-300)
p.qgammaSml(1e-300,1e-50, a2= seq(100,360, by=4), a3=seq(350,1500, by=10))
showProc.time()
## The "upper" problematic corner:
p.qgammaSml(1e-19, 1e-3, a2=NULL,a3=NULL, ylim=c(.1,20))
p.qgammaSml(1e-19, 1e-3, a2=seq(1,12, by=.04), ylim=c(.1,20),a3=NULL,k.lab=10)
## now shows the problem (quite well):
## could it be in pgamma()'s inaccuracy, leading to qgamma() bias ?
aa <- c(seq(9, 22, by=0.25),seq(22.3,40,by=0.4))
caa <- formatC(range(aa))
sfsmisc::mult.fig(2)
curve(pgamma(x, sh=aa[1]), 0.5, 20, log = 'xy', ylim = c(1e-60, .2),
main = sprintf("pgamma(x, a) for a in [%s,%s]", caa[1],caa[2]))
for(sh in aa) curve(pgamma(x, sh), add = TRUE, col=2)
abline(h=c(1e-15), col="light blue", lty=2)
curve(pgamma(x, sh=aa[1]), 0.5, 20, log = 'xy', ylim = c(1e-15, .8),
main = sprintf("pgamma(x, a) for a in [%s,%s]", caa[1],caa[2]))
for(sh in aa) curve(pgamma(x, sh), add = TRUE, col=2)
## the "border curve" between "Pearson" and "Continued fraction (upper tail)"
## in pgamma.c :
curve(pgamma(max(1,x), x), add = TRUE, col=4)
## ==> pgamma() is perfect here {series expansion up to eps_C accuracy}!
aa <- c(seq(9, 22, by=0.25),seq(22.3,40.4,by=0.4))
p.qgammaSml(1e-24, 1e-5, a1=aa, a2=NULL,a3=NULL, ylim=c(.8,8))
## -------- save the above?
aa1 <- c(aa,seq(40.5,90, by=0.5))
p.qgammaSml(1e-60, 1e-5, a1=aa1, a2=NULL,a3=NULL, ylim=c(.9, 16))
aa2 <- c(aa1, seq(91,150, by= 1))
p.qgammaSml(1e-90, 1e-5, a1=aa2, a2=NULL,a3=NULL, ylim=c(.9, 35))
aa3 <- c(aa2, seq(150,250, by= 2), seq(253, 400, by=5))
p.qgammaSml(1e-200, 1e-5, a1=aa3, a2=NULL,a3=NULL, ylim=c(.9, 100))
p.qgammaSml(1e-200, 1e-5, a1=aa3, a2=NULL,a3=NULL, ylim=c(.9, 200),k.lab=9e9)
p.qgammaSml(1e-60, 1e-5, a1=aa3, a2=NULL,a3=NULL, ylim=c(.9, 200),k.lab=9e9)
showProc.time()
## lower a \> 10
curve(qgamma(x, 19), 1e-14, 1e-9, log='x')
curve(qgamma(x, 18), 1e-14, 1e-9, log='x')
curve(qgamma(x, 15), 1e-11, 5e-9, log='x')
curve(qgamma(x, 13), 5e-10, 1e-8, log='x')
curve(qgamma(x, 11), 1e-8, 5e-8, log='x')
curve(qgamma(x, 10.5), 4.2e-8, 6e-8, log='x')
curve(qgamma(x, 10.3), 6e-8, 7e-8, log='x')
curve(qgamma(x, 10.2), 7.1e-8, 7.6e-8, log='x')
curve(qgamma(x, 10.15),7.7e-8, 7.9e-8, log='x')
curve(qgamma(x, 10.14),7.88e-8,7.92e-8, log='x',n=10001)
## no more problems for smaller a!! here:
curve(qgamma(x, 10.13), 1e-10, 5e-4, log='x',n=20001)
curve(qgamma(x, 10.12), 1e-10, 5e-4, log='x',n=20001)
curve(qgamma(x, 10.1), 1e-10, 5e-4, log='x',n=20001)
showProc.time()
##--- the "+Inf" / premature "0" case:
curve(qgamma(x, 155, log=TRUE), -1500, 0, log='y', n=2001,col=2)
curve(qgamma(x, 1e3, log=TRUE), -1500, 0, log='y', n=2001,col=2)
## now works, but slowly and with kink
curve(qgamma (x, 1e5, log=TRUE), -3e5, 0, log='y', n=2001,col=2,lwd=3)
curve(qgammaAppr(x, 1e5, log=TRUE), add = TRUE, n=2001, col="blue",lwd=.4)
## --- curves are almost "identical"
## ===> the kink *does* come from the initial approx... hmm
## still "identical"
curve(qgamma (x, 1e4, log=TRUE), -3e4, 0, log='y', n=2001,col=2)
curve(qgammaAppr(x, 1e4, log=TRUE), add = TRUE, n=2001, col="tomato3")
## now see some difference (approx. has kink at ~ -165)
curve(qgamma (x, 100, log=TRUE), -200, 0, log='y', n=2001,col=2)
curve(qgammaAppr(x, 100, log=TRUE), add = TRUE, n=2001, col="tomato3")
##
(kk <- 100 * 2/1.24)# 161.29
curve(qgamma (x, 100, log=TRUE), -1.1*kk, -.95*kk, log='y', n=2001,col=2)
curve(qgammaAppr(x, 100, log=TRUE), add = TRUE, n=2001, col="tomato3")
abline(v = -kk, col='blue', lty=2)# exactly: kink is at a * 2 / 1.24 = a / .62
curve(qgammaAppr(x - 100/.62, 100,log=TRUE), -1e-3, +1e-3)
showProc.time()
p.qgammaLog <- function(alpha, xl.f = 1.5, xr.f = 0.4, n = 2001)
{
## Purpose:
## ----------------------------------------------------------------------
## Arguments:
## ----------------------------------------------------------------------
## Author: Martin Maechler, Date: 30 Mar 2004, 18:44
kk <- -alpha / .62 # = (alpha * 2) / (-1.24)
curve(qgamma(x, alpha, log=TRUE), xl.f*kk, xr.f*kk, log='y',
n=n, col=2, lwd=3.6, lty = 4,
main= paste("qgamma(x, alpha=",formatC(alpha,digits=10),", log = TRUE)"))
lines(kk, qgamma(kk, alpha, log=TRUE), type = 'h', lty = 3)
curve(qgamma (exp(x), alpha), add = TRUE, col="orange", n=n, lwd= 2)
curve(qgammaAppr(x, alpha, log=TRUE), add = TRUE, col=3, n=n,lwd = .4)
}
showProc.time()
p.qgammaLog(25)
p.qgammaLog(16)# ~ [-25, -20]
p.qgammaLog(12, 1.2, 0.8)# small problem remaining
p.qgammaLog(11, 1.2, 0.8)# even smaller
p.qgammaLog(10.5, 1.1, 0.9)# even smaller
p.qgammaLog(10.25, 1.1, 0.9)# even smaller
## 2019-08: __nothing__ visible from here on:
p.qgammaLog(10.18, 1.02, 0.98)# even smaller
p.qgammaLog(10.15, 1.02, 0.98)# even smaller
p.qgammaLog(10.14, 1.001, 0.999)# even smaller
p.qgammaLog(10.139, 1.0002, 0.9998)#
p.qgammaLog(10.138, 1.0002, 0.9998)#
p.qgammaLog(10.137, 1.00001, 0.99999)#
p.qgammaLog(10.13699, 1.0000001, 0.9999999)#
p.qgammaLog(10.1369899, 1.0000001, 0.9999999)#
p.qgammaLog(10.1369894, 1.0000001, 0.9999999)#
p.qgammaLog(10.1369893, 1.0000001, 0.9999999)# even smaller at -16.34998
showProc.time()
##-- here is the boundary --- for 64-bit AMD Opteron ---
## and for 32-bit AMD Athlon
p.qgammaLog(10.1369892, 1.0000001, 0.9999999)# no more
p.qgammaLog(10.136989, 1.0000001, 0.9999999)#
p.qgammaLog(10.136988, 1.0000001, 0.9999999)#
p.qgammaLog(10.136985, 1.0000001, 0.9999999)#
p.qgammaLog(10.13698, 1.0000001, 0.9999999)#
p.qgammaLog(10.13697, 1.0000001, 0.9999999)#
p.qgammaLog(10.13695, 1.0000001, 0.9999999)#
p.qgammaLog(10.1368, 1.000001, 0.999999)#
p.qgammaLog(10.1365, 1.000001, 0.999999)#
p.qgammaLog(10.136, 1.000001, 0.999999)#
p.qgammaLog(10.125, 1.1, 0.9)# --- see it now
p.qgammaLog(10, 1.2, 0.8)
p.qgammaLog(9)
showProc.time()
## For large alpha: show difference to see problem better
## ---> for alpha >= 10, the x problem starts *roughly* at x = -0.8*alpha
##
sfsmisc::mult.fig(2)
curve(qgammaAppr(x, 5, log=TRUE), - 8.1, -8, n=2001)
curve(qgammaAppr(x- 5/.62, 5, log=TRUE), -1e-15, 0)
## is the kink from pgamma() ? : no: this looks fine,
curve(pgamma(x, 1e5, log=TRUE), 1, 2e5, log='x', n=2001,col=2)
## and this does too:
curve( dgamma(x, 1e5), .5e5, 2e5); par(new=TRUE)
curve( dgamma(x, 1e5, log=TRUE), .5e5, 2e5, col=2, yaxt="n")
axis(4,col.axis=2); par(new=TRUE)
curve( pgamma(x, 1e5), .5e5, 2e5, n=2001, col=3); par(new=TRUE)
curve( pgamma(x, 1e5, log=TRUE), .5e5, 2e5, n=2001, col=4); par(new=TRUE)
curve(-pgamma(x, 1e5, log=TRUE,lower=FALSE), .5e5, 2e5, n=2001, col=4)
## all looking nice
x <- 10^seq(2,6, length=4001)
qx <- qgamma(pgamma(x, 1e5, log=TRUE), 1e5, log=TRUE)
plot(x, qx, type ='l', col=2, asp = 1); abline(0,1, lty=3)
showProc.time()
###------------- Approximations of qgamma() ------
##
## source("/u/maechler/R/MM/NUMERICS/dpq-functions/qchisqAppr.R")
##--> qchisqAppr()
##--> qchisqWH [ = Wilson Hilferty ]
##--> qchisqKG [ = Kennedy & Gentle's improvements "a la AS 91" ]
## dyn.load("/u/maechler/R/MM/NUMERICS/dpq-functions/qchisq_appr.so")
## Consider the two different implementations of
## lgamma1p(a) := lgamma(1+a) == log(gamma(1+a) == log(a*gamma(a)) "stable":
lseq <- if(requireNamespace("sfsmisc")) sfsmisc::lseq else
function(from, to, length) exp(seq(log(from), log(to), length.out = length))
if(require("Rmpfr")) { ##---------------- MPFR numbers -------------------------
.mpfr.all.eq <- Rmpfr::all.equal
AllEq <- function(target, current, ...)
.mpfr.all.eq(target, current, ...,
formatFUN = function(x, ...) Rmpfr::format(x, digits = 9))
print(gammaE <- Const("gamma",200)); pi. <- Const("pi",200)
print(a0 <- (gammaE^2 + pi.^2/6)/2)
print(psi2.1 <- -2*zeta(mpfr(3,200)))# == psigamma(1,2) =~ -2.4041138
print(a1 <- (psi2.1 - gammaE*(pi.^2/2 + gammaE^2))/6)
x <- lseq(1e-30, 0.8, length = if(doExtras) 1000 else 125)
x. <- mpfr(x, 200)
xct. <- log(x. * gamma(x.)) ## using MPFR arithmetic .. no overflow ...
xc2. <- log(x.) + lgamma(x.)## (ditto)
print(AllEq(xct., xc2., tol = 0)) # 3.15779......e-57
xct <- as.numeric(xct.)
stopifnot(exprs = {
AllEq(xct., xc2., tol = 1e-45)
AllEq(xct , xc2., tol = 1e-15)
##
all.equal(lgamma1p(x), lgamma1p(x, tol= 1e-16), tol=0)
## -> no difference; i.e., default tol = 1e-14 seems fine enough!
})
showProc.time()
m.appr <- cbind(log(x*gamma(x)), lgamma(1+x), log(x) + lgamma(x),
lgamma1p.(x, k=1, cut=3e-6),
lgamma1p.(x, k=2, cut=1e-4),
lgamma1p.(x, k=3, cut=8e-4),
lgamma1p(x))#, tol= 1e-14), # = default
eMat <- m.appr - xct # absolute error
## Relative errors:
str(reMat. <- m.appr/xct. - 1)
str(reMat <- as(reMat., "array")) # as(., "matrix") fails in older versions
matplot(x, eMat , log="x", type="l", lty=1) #-> problematic log(x) + lgamma(x) for "large"
matplot(x, abs( eMat), log="xy", type="l", lty=1) #-> but good for small; lgamma1p is much better
matplot(x, abs(reMat), log="xy", type="l", lty=1)
abline(v= 3.47548562941137e-08, col = "gray80", lwd=3)#<- the cutoff value of lgamma1p()
##---> should use earlier cutoff!
## zoom in:
matplot(x, abs(reMat), log="xy", type="l", col=1:7, lty=1,
lwd=2, xlim=c(8e-9, 1e-3), ylim = c(1e-18, 1e-7), axes=FALSE, frame=TRUE,
main = expression(lgamma1p(x) == log(Gamma(x+1)) ~~~ "approximations"
~~~ abs(rel.Err(.))))
eaxis(1); eaxis(2)
abline(v= 3.47548562941137e-08, col = "gray80", lwd=3)#<- the cutoff value of lgamma1p()
abline(h= c(1,2,4)*.Machine$double.eps, lty=3, col="skyblue")
legend("topright", col=1:7, lty=1,lwd=2,
c("log(x*gamma(x))", "lgamma(1+x)", "log(x) + lgamma(x)",
"lgamma1p.(x, k=1, c=3e-6)",
"lgamma1p.(x, k=2, c=1e-4)",
"lgamma1p.(x, k=3, c=8e-4)",
"lgamma1p(x)"), bty="n", ncol=2)
abline(v = c(3e-6, 1e-4, 8e-4), col=4:6, lty=2, lwd=1/2)
## FIXME: do the same for the lgaamma1p_series()
## rm(x., xct., xc2., reMat., eMat, AllEq)
detach("package:Rmpfr")
showProc.time()
} ## if( MPFR ) ----------------------------------------------------------------
## ../R/qchisqAppr.R -- talks about the "small shape" qgamma() approximation
## ----------------- --> .qgammaApprBnd() :
curve(.qgammaApprBnd, 1e-18, 1e-15, col=2)
abline(h=0, col="gray70", lty=2)
eps.c <- .Machine$double.eps
axis(3,at=(1:3)* eps.c,
label=expression(epsilon[c], 2*epsilon[c], 3*epsilon[c]))
(rt.b <- uniroot(.qgammaApprBnd, c(1,3)*eps.c, tol=1e-12))
rt.b$root ## 3.954775e-16
rt.b$root / eps.c ## 1.781072
##==> for a < 1.781*eps, bnd > 0 ==> we have log(p) < bnd for all p
## otherwise, we should effectively 'test'
curve(.qgammaApprBnd, 1e-16, 1e-10, log="x", col=2)
showProc.time()
## source ("/u/maechler/R/MM/NUMERICS/dpq-functions/beta-gamma-etc/qgamma-fn.R")
## ##--> qchisqAppr.R() -- which has 'kind = ' argument!
## ##--> qgamma.R()
p.qchi.appr <-
function(x, qm= { m <- cbind(qchisq(x, df, log=TRUE),
sapply(knds, function(kind)
qchisqAppr.R(x,df,log=TRUE,kind=kind)))
colnames(m) <- c("True", "default", knds[-1])
m },
df,
knds = list(NULL,"chi.small", "WH", "p1WH", "df.small"),
call = match.call(), main = deparse(call), log = "y", ...)
{
## Purpose:
## ----------------------------------------------------------------------
## Arguments:
## ----------------------------------------------------------------------
## Author: Martin Maechler, Date: 25 Mar 2004, 22:08
col <- c(2,1,3:6)
lty <- c(1,3,1,1,1,1)
lwd <- c(2,2,1,1,1,1)
matplot(x, qm, col=col, lty=lty, lwd=lwd, log = log,
type = 'l', main = main, ...)
y0 <- c( .02, .98) %*% par('usr')[3:4]
if(par("ylog")) y0 <- 10^y0
legend(min(x), y0, colnames(qm), col=col, lty=lty, lwd=lwd)
invisible(list(x=x, qm=qm, call=match.call()))
}
pD.chi.appr <- function(pqr, err.kind=c("relative", "absolute"),
type = "l", log = "y",
lwds = c(2, rep(1, k-1)),
cols = seq(along=lwds),
ltys = rep(1,k),
...)
{
## Purpose: Plot Difference from "True" qchisq()
## ----------------------------------------------------------------------
## Arguments: pqr: a list as resulting from p.chi.appr()
## ----------------------------------------------------------------------
## Author: Martin Maechler, Date: 31 Mar 2004, 09:38
err.kind <- match.arg(err.kind)
if(!is.list(pqr) || !is.numeric(k <- ncol(pqr$qm)-1) || k <= 0)
stop("invalid first argument 'pqr'")
with(pqr, {
err <- abs(if(err.kind == "relative")
(1- qm[,-1] / qm[,1]) else (qm[,-1] - qm[,1]))
matplot(x, err, ylab = "",
main = paste(err.kind,"Error from", deparse(call)),
type= type, log= log, lty=ltys, col=cols, lwd=lwds, ...)
legend(par("xaxp")[1], par("yaxp")[2], colnames(qm)[-1],
lty=ltys, col=cols, lwd=lwds)
})
}
## if(FALSE) # you can manually set
## do.pqchi <- TRUE # before source()ing this file
## if(!exists("do.pqchi") || !is.logical(do.pqchi))
## do.pqchi <- interactive()
## if(do.pqchi) { #------------- FIXME look at speed up or cache indeed !?? <<<<<
## pqFile <- "/u/maechler/R/MM/NUMERICS/dpq-functions/pqchi.rda"
## ## ls -l ................... 1325446 Nov 2 2009 pqchi.rda
## if(file.exists(pqFile)) {
## attach(pqFile) ## it loads more than we create here __FIXME__
## print(ls(2, all.names=TRUE))
## ## [1] "pq.1" "pq.25" "pq.25." "pq.31" "pq.33" "pq.33." "pq.33.2"
## ## [8] "pq.33.3" "pq.33.4" "pq.5" "pq.5." "pq1" "pq1." "pq1.95"
## ## [15] "pq1.95." "pq1.95.2" "pq10" "pq10." "pq10.2" "pq100" "pq2"
## ## [22] "pq2." "pq2.05" "pq2.05." "pq2.05.2" "pq2.5" "pq2.5." "pq2.5.2"
## ## [29] "pq200" "pq2L" "pq4" "pq4." "pq4.2" "pqFile"
## }
showProc.time()
x <- seq(-300, -.01, length=501)
all.equal(qchisqAppr.R(x, 200, log=TRUE),
qchisqAppr (x, 200, log=TRUE),tol=0)
## 4.48 e-16 / TRUE (Opteron)
all.equal(qchisqAppr.R(x, 2, log=TRUE),
qchisqAppr (x, 2, log=TRUE),tol=0)
## 3.90 e-16 / TRUE (Opteron)
all.equal(qchisqAppr.R(x, 0.1, log=TRUE),
qchisqAppr (x, 0.1, log=TRUE),tol=0)
## 7.15 e-15 / 1.179e-8 !!!!! (Opteron)
pq200 <- p.qchi.appr(x = seq(-300, -.01, length=501), df = 200)
pq100 <- p.qchi.appr(x = seq(-160, -.01, length=501), df = 100)
## after (slow) computing, quickly repeat:
with(pq200, p.qchi.appr(x=x, qm=qm, call=call))
with(pq100, p.qchi.appr(x=x, qm=qm, call=call))
## this "hangs forever" -- before I introduced 'maxit' (for 'nu.small'):
pq10 <- p.qchi.appr(x = seq(-12, -.005, length=501), df = 10)
## want to see the jump:
pq10. <- p.qchi.appr(x = seq(-10, -4, length=501), df = 10)
pq10.2<- p.qchi.appr(x = seq(-8.5,-7.5, length=501), df = 10)
with(pq10.2, p.qchi.appr(x=x, qm=qm, call=call))
pq2.5 <- p.qchi.appr(x = seq(-3.4, -.01, length=501), df = 2.5)
pq2.5. <- p.qchi.appr(x = seq(-2.1, -1.8, length=901), df = 2.5)#the jump
## what about p1WH (which is fantastic for df=2)?
pq2.5.2<- p.qchi.appr(x = seq(-0.5, -1e-3, length=901), df = 2.5)
with(pq2.5, p.qchi.appr(x=x, qm=qm, call=call))
with(pq2.5.2, p.qchi.appr(x=x, qm=qm, call=call))
pD.chi.appr(pq2.5.2)# nothing special
pq2.05 <- p.qchi.appr(x = seq(-3.4, -.01, length=501), df = 2.05)
pq2.05. <- p.qchi.appr(x = seq(-2.5, -1.5, length=901), df = 2.05)#the jump
## ^^ the jump from chi.small to WH is much too late here
## what about p1WH (which is fantastic for df=2)?
pq2.05.2<- p.qchi.appr(x = seq(-0.4, -1e-5, length=901), df = 2.05)
pD.chi.appr(pq2.05.2) # p1WH is starting to become better
# and the jump (WH -> p1WH) is too late
with(pq2.05, p.qchi.appr(x=x, qm=qm, call=call))
with(pq2.05.2, p.qchi.appr(x=x, qm=qm, call=call))
## Here, all are 'ok' (but "nu.small"):
pq2L <- p.qchi.appr(seq(-300, -.01, length=201), df = 2)
pq2 <- p.qchi.appr(x = seq(-5, -.001, length=501), df = 2)
pq2. <- p.qchi.appr(x = seq(-2.5, -1, length=901), df = 2)
with(pq2., p.qchi.appr(x=x, qm=qm, call=call))
pq4 <- p.qchi.appr(x = seq(-8, -0.01, length = 501), df = 4)
summary(warnings())
pq4.2 <- p.qchi.appr(x = seq(-0.1, -1e-04, length = 901), df = 4)
pq1.95 <- p.qchi.appr(x = seq(-3., -.01, length=501), df = 1.95)
pq1.95. <- p.qchi.appr(x = seq(-2.2, -1.5, length=901), df = 1.95)#the jump -1.57
## ^^ the jump from chi.small to WH is *much* too late here
## what about p1WH (which is fantastic for df=2)?
pq1.95.2<- p.qchi.appr(x = seq(-0.4, -1e-7, length=901), df = 1.95)
pD.chi.appr(pq1.95.2) # p1WH is starting to become better
# and the jump (WH -> p1WH) is too late
with(pq1.95, p.qchi.appr(x=x, qm=qm, call=call))
with(pq1.95.2, p.qchi.appr(x=x, qm=qm, call=call))
pq1 <- p.qchi.appr(x = seq(-4, -.001, length=501), df = 1)
pq1. <- p.qchi.appr(x = seq(-1.8, -.5, length=901), df = 1)
with(pq1., p.qchi.appr(x=x, qm=qm, call=call))
pq.5 <- p.qchi.appr(x = seq(-1.5, -.001, length=501), df = 0.5)
pq.5. <- p.qchi.appr(x = seq(-0.8, -0.2, length=901), df = 0.5, ylim=c(.04,.6))
pq.33 <- p.qchi.appr(x= seq(-0.9, -.001,length=501), df= 0.33)
pq.33. <- p.qchi.appr(x= seq(-0.4, -0.02,length=901), df= 0.33)
pq.33.2<- p.qchi.appr(x= seq(-0.4, -0.2, length=901), df= 0.33, ylim=c(.15,.60))
with(pq.33.2, p.qchi.appr(x=x, qm=qm, call=call,ylim=c(.15,.45)))
pq.33.3<- p.qchi.appr(x= seq(-0.4, -0.005, length=901), df= 0.33, ylim=c(.15, 4.00))
with(pq.33.3, p.qchi.appr(x=x, qm=qm, call=call))#,ylim=c(.15, 8)))
pq.33.4<- p.qchi.appr(x= seq(-0.003, -1e-6, length=901), df= 0.33,ylim=c(5,25))
with(pq.33.4, p.qchi.appr(x=x, qm=qm, call=call,ylim=c(5,25)))
## nu <= 0.32 is the "magic" border of Best & Roberts
pq.31 <- p.qchi.appr(x = seq(-0.45,-.010, length=501), df = 0.31)
with(pq.31, p.qchi.appr(x=x, qm=qm, call=call))
pq.25 <- p.qchi.appr(x = seq(-0.3, -0.02, length=901), df = 0.25)
pq.1 <- p.qchi.appr(x = seq(-0.16, -0.01, length=901), df = 0.1)
with(pq.1, p.qchi.appr(x=x, qm=qm, call=call))
showProc.time()
## if(!file.exists(pqFile)) # don't overwrite for now (as it contains pq2L ,
## save(list=ls(pat="^pq"), file = pqFile)
##}## end if(do.pqchi){ only if interactive } ======================================
pD.chi.appr(pq2L, "abs")
pD.chi.appr(pq2L, "rel")
## --> want only much smaller x-range:
pD.chi.appr(pq2,"abs")#--> fantastic p1WH
pD.chi.appr(pq2) # (ditto)
pD.chi.appr(pq4.2)# p1WH: only at very right
pD.chi.appr(pq4.2, xlim=c(-.016,0))# p1WH: only at very right
## no Newton step here, eg:
(qgA100 <- qgammaAppr(1e-100, 100))
(qg.100 <- qgamma (1e-100, 100))
all.equal(qgA100, qg.100)
## too much different
dgamma(1e-100, 100, log = TRUE)# -23154.7 i.e., "non-log" is 0
qgamma.R(1e-100, 100, verbose = TRUE)#-> final Newton fails!
## but here, the final Newton iterations do work :
x <- qgamma.R(log(1e-100), 100, log = TRUE, verbose = TRUE)
pgamma(x, 100) # = 1e-100 ! perfect
showProc.time()
###--> Use this to devise an improved final Newton algorithm !!!
## From: Prof Brian Ripley <ripley@stats.ox.ac.uk>
## To: skylab.gupta@gmail.com
## Cc: R-bugs@biostat.ku.dk, r-devel@stat.math.ethz.ch
## Subject: Re: [Rd] qgamma inaccuracy (PR#12324)
## Date: Tue, 12 Aug 2008 20:50:50 +0100 (BST)
## This is a really extreme usage. AFAICS the code works well enough down to
## shape=1e-10 or so, e.g.
qgamma(1e-10, 5e-11, lower.tail=FALSE)
## [1] 0.08237203
## in R 2.9.. .. 2.10.0 -- with an accuracy warning {which is *wrong*!}
## I would be interested to know what substantive problem you were trying to
## solve here that required such values.
## I am pretty sure that a completely different algorithm will be required.
## MM: It looks like this is basis for a new algo:
a <- 1e-14
gammaE <- 0.57721566490153286060651209008240243079 # Euler's gamma
curve(pgamma(x, a, lower.tail=FALSE)/a + log(x) + gammaE, 1e-300, 1e-1, log="",n=1000)
curve(pgamma(x, a, lower.tail=FALSE)/a + log(x) + gammaE - x, 1e-300, 1e-1, log="",n=1000)
## ==> Q = 1 - P = pgamma(x,a, lower=FALSE) ~= a*(-log(x) - gammaE + x - x^2/4)
## i.e., Q ~= -a(log(x) + gammaE { -x + x^2/4 }
## -Q/a - gammaE ~= log(x) { -x + x^2/4 }
## ==> x ~= exp(- (Q/a + gammaE))
## e.g., example below:
Q <- 1e-100; a <- 5e-101
## MM: Find inverse :
str(r.a <- uniroot(function(x) pgamma(x,a,lower.tail=FALSE) - Q,
int = c(0.01, 0.1), tol=1e-20))
dput(x0 <- r.a$root) ## 0.0823720296207203
(x1 <- exp(- (Q/a + gammaE)))## 0.07598528 .. not so good
qgammaApprSmallP(Q, a, lower.tail=FALSE)## ~= 0.07598528 -- the same!!
pgamma(x0, a, lower.tail=FALSE) ## 1.00000e-100.
pgamma(x1, a, lower.tail=FALSE) ## 1.03728e-100 ... hmm "close"
##
## MM: -- now look at the bigger picture
p.qg.2a <- function(l2x.min= -15, l2x.max = -100, n = 501,
do.offset = FALSE,
type = "o", log = "x", cex = 0.6, ...) {
x.log <- any("x" == strsplit(log,"")[[1]])
x <- if(x.log) 2^seq(l2x.min, l2x.max, length=n)
else seq(2^l2x.min,2^l2x.max, length=n)
if(do.offset)
plot(x, qgamma(2*x, x, lower.tail=FALSE) - 0.0823720296206873,
type=type, cex=cex, log=log, ...)
else plot(x, qgamma(2*x, x, lower.tail=FALSE),
type=type, cex=cex, log=log, ...)
}
p.qg.2a() # was "very bad" in R <= 2.10.0, now --> 0.082372... "perfect" smooth
## still a little ---but acceptable--- remaining inaccuracy ...zooming in:
p.qg.2a(-43,-55, do.offset=TRUE)
p.qg.2a(,-1024)
p.qg.2a(,-1024, log="", pch=".")## linear in x !!
## zoom in at the limit
p.qg.2a(-30,-1024, do.offset=TRUE, ylim = 1e-11*c(-1,1))
p.qg.2a(-33,-1024, do.offset=TRUE, ylim = 1e-12*c(-1,1))
p.qg.2a(-33,-1024, do.offset=TRUE, ylim = 1e-13*c(-1,1))
a <- 2^-(7:900)
qg <- qgamma(2*a, a, lower.tail=FALSE)
re <- 1-pgamma(qg, a, lower.tail=FALSE)/(2*a)
plot(a, re, log="x", type="b", col=2)
stopifnot(abs(re) < 2e-12) # but really, *should be a bit better
showProc.time()
## For completeness we may write that in due course, but for now (R 2.7.2) I
## suggest just issuing a warning for miniscule 'shape'.
## On Thu, 7 Aug 2008, skylab.gupta@gmail.com wrote:
## > Full_Name:
## > Version: 2.7.1 (2008-06-23)
## > OS: windows vista
## > Submission from: (NULL) (216.82.144.137)
## >
## >
## > Hello,
## >
## > I have been working with various probability distributions in R, and it seems
## > the gamma distribution is inaccurate for some inputs.
## > For example, qgamma(1e-100, 5e-101, lower.tail=FALSE) gives: 1.0. However, it
## > seems this is incorrect; I think the correct answer should be
## > 0.082372029620717283. When I check these numbers using pgamma, I get:
(qg <- qgamma(1e-100, 5e-101, lower.tail=FALSE))# 1 (wrong, originally)
## 0.08237203 now (2009-11-04), i.e. ok
pgamma(qg, 5e-101, lower.tail=FALSE)# now -> 1e-100 : ok
pgamma(0.082372029620717283,5e-101, lower.tail=FALSE)
## 1.0000000000000166e-100.
RE.pqgamma <- function(p, shape, lower.tail = TRUE, log.p = FALSE) {
## Relative Error of pgamma(qgamma(*), ..):
1 - pgamma(qgamma(p, shape, lower.tail=lower.tail, log.p=log.p),
shape=shape, lower.tail=lower.tail, log.p=log.p) / p
}
RE.qpgamma <- function(q, shape, lower.tail = TRUE, log.p = FALSE) {
## Relative Error of qgamma(pgamma(*), ..):
1 - qgamma(pgamma(q, shape, lower.tail=lower.tail, log.p=log.p),
shape=shape, lower.tail=lower.tail, log.p=log.p) / q
}
## Ok, how extreme can we get -- let a := alpha := shape --> 0 :
x <- 1e-100
a <- 2^-(7:300)# is still "ok":
plot(a, (re <- RE.pqgamma(x, a, lower.tail=FALSE)), type="b", col=2, log="x")# oops!
a <- 2^-(7:400)
plot(a, (re <- RE.pqgamma(x, a, lower.tail=FALSE)), type="b", col=2, log="x")# oops!
## Oops!
## but, it is clear
qgamma(x, 2^-400, lower.tail=FALSE)## is exactly 0
## -> it goes to 0 quickly . .. zooming in:
curve(qgamma(1e-100, x, lower.tail=FALSE), 1e-110, 1e-70, log="xy", col=2, axes=FALSE)
eaxis(1);eaxis(2)
## from when on is it exactly 0:
uniroot(function(u) qgamma(1e-100, 2^u,lower.tail=FALSE)-1e-315, c(-400, -300))$root
## -341.6941
## => use
a <- 2^-(7:341)
plot(a, (re <- RE.pqgamma(x, a, lower.tail=FALSE)),
type="l", col=2, log="x")# small glips
## zoom in:
curve(abs(RE.pqgamma(1e-100, x, lower.tail=FALSE)), 2^-341, 1e-90, log="xy", n=2000)
curve(RE.pqgamma(1e-100, x, lower.tail=FALSE), 1e-100, 10e-100, n=2000)
curve(RE.pqgamma(1e-100, x, lower.tail=FALSE), 4e-100, 6e-100, n=2000)
## Ok: at least here is a problem
RE.pqgamma(1e-100, 5e-100, lower.tail=FALSE)# -0.1538
## more general
curve(RE.pqgamma(x, 5*x, lower.tail=FALSE), 1e-100, 1e-20, n=10000, log="x")
## problem *everywhere* , starting quite early: (lesser problem at ~ 1e-16 !)
curve(RE.pqgamma(x, 5*x, lower.tail=FALSE), 1e-25, 1e-15, n=1000, log="x")
curve(RE.pqgamma(x, 5*x, lower.tail=FALSE), 1e-21, 10e-21,n=1000)
curve(RE.pqgamma(x, 5*x, lower.tail=FALSE), 2e-21, 6e-21, n=1000)
curve(RE.pqgamma(x, 5*x, lower.tail=FALSE), 4e-21, 4.5e-21, n=1000)
## and indeed, it's qgamma() that jumps here
curve(qgamma(x, 5*x, lower.tail=FALSE), 4e-21, 4.5e-21, ylim=c(.97, 1.1))
## well, looking at pgamma(), finally reveals the buglet is there first:
## There's a jump at x = 1 !!
curve(pgamma(x, 1e-30, lower=FALSE), .9999, 1.0001)
curve(pgamma(x, 1e-20, lower=FALSE), .9999, 1.0001)
curve(pgamma(x, 1e-17, lower=FALSE), .9999, 1.0001)
curve(pgamma(x, 1e-15, lower=FALSE), .9999, 1.0001)
curve(pgamma(x, 1e-13, lower=FALSE), .9999, 1.0001)
curve(pgamma(x, 1e-12, lower=FALSE), .9999, 1.0001)# still clearly visible
curve(pgamma(x, 1e-11, lower=FALSE), .9999, 1.0001)# barely visible
## for larger alpha == shape, must zoom in more and more:
curve(pgamma(x, 1e-11, lower=FALSE), .99999, 1.00001)#
curve(pgamma(x, 1e-10, lower=FALSE), .999999, 1.000001)
curve(pgamma(x, 1e-9, lower=FALSE), .9999999, 1.0000001)
curve(pgamma(x, 1e-8, lower=FALSE), .99999999, 1.00000001)
curve(pgamma(x, 1e-7, lower=FALSE), .999999999, 1.000000001)
curve(pgamma(x, 1e-6, lower=FALSE), .9999999999, 1.0000000001)
curve(pgamma(x, 1e-5, lower=FALSE), .99999999999, 1.00000000001)
curve(pgamma(x, 1e-4, lower=FALSE), .999999999999, 1.000000000001)
## now we get close to noise level:
curve(pgamma(x, 1e-3, lower=FALSE), .9999999999999, 1.0000000000001)
curve(pgamma(x, 1e-2, lower=FALSE), .99999999999999, 1.00000000000001)
showProc.time()
del.pgamma <- function(a, eps = 1e-13)
{
## Purpose: *relative* jump size at x = 1 of pgamma(x, a, lower=FALSE)
## ----------------------------------------------------------------------
## Author: Martin Maechler, Date: 5 Nov 2009, 16:08
stopifnot(a > 0, length(a) == 1, eps > 0, length(eps) == 1)
pp <- pgamma(1+c(-eps,eps), a, lower.tail = FALSE)
(pp[2] - pp[1])*2/(pp[2] + pp[1])
}
a <- lseq(1e-300, 1e-3, length=400)
dpa <- sapply(a, del.pgamma)
plot(a, -dpa, log="xy", type="l", col=2, yaxt="n");eaxis(2)
## ok, it remains constant all the way to 1e-300
## --> focus
a <- lseq(1e-40, 1e-5, length=400)
dpa <- sapply(a, del.pgamma)
plot(a, -dpa, log="xy", type="l", col=2, axes=FALSE)
eaxis(1, at = 10^-seq(5,40, by=5))
eaxis(2)
xm <- .Machine$double.xmin
pgamma(xm, shape= 1e-20)# is "practically 1" --> *most* qgamma() will be exactly 0
## how close to 1 ? ---> use upper_tail [possibly on log scale:]
pgamma(xm, shape= 1e-20, lower.tail=FALSE) # 7.078192e-18
pgamma(xm, shape= 1e-20, lower.tail=FALSE, log=TRUE) # -39.48951
## Where is the 'boundary' (from which on qgamma() must return 0, since it can't give
## xm = 2.2e-308 ) ?
a <- 2^-(7:1030)
plot(a, (p <- pgamma(xm, a, lower.tail=FALSE, log=TRUE)),
cex=.5,type ="l", col=2, log="x")
summary(lm. <- lm(p ~ log(a) + a + I(a^2))) ## coeff. of log(a) *is* 1
summary(lm. <- lm(p ~ offset(log(a)) + a + I(a^2)))
p.l <- pgamma(xm, a, lower.tail=FALSE, log=TRUE) - log(a)
dput(mean(tail(p.l))) ## 6.56218869790132
##=> pgamma(xm, a) ~= log(a) + 6.5621886979 +
## ok, to get this better, now need different a:
al <- seq(1e-300, 1e-3, length=200)
plot(al, (pl <- pgamma(xm, al, lower.tail=FALSE, log=TRUE)) - (log(al)+6.56218869790132),
cex=.5,type ="l", col=2)
summary(lm.. <- lm(pl ~ offset(log(al) + 6.56218869790132) + 0 + al + I(al^2)))
confint(lm..)
coef(lm..)
## al I(al^2)
## -353.8587 20750.8205
##=> pgamma(xm, a) ~= log(a) + 6.5621886979 - 353.858745 * a
## > Similarly, for example: -- now (2009-11-04) ok
qgamma(1e-100,0.005,lower.tail=FALSE) # = 109.36757177 instead of 219.5937..
pgamma(109.36757177007101, 0.005, lower.tail=FALSE)# = 1.4787306506694758e-52.
## > This looks completely wrong. The correct value, I think, should be
## > 219.59373661415756. In fact,
pgamma(219.59373661415756, 0.005, lower.tail=FALSE)# = 9.9999999999999558e-101.
## >
## > In fact, when I do the following in R, the results are completely wrong,
## >
a <- 5*10^-(1:40)
z1 <- qgamma(1e-100,a,lower.tail=FALSE)
(y <- pgamma(z1,a,lower.tail=FALSE))
## The value of y that I get should be close to 1e-100, but they are not:
## [1] 1.000000e-100 1.871683e-51 1.478731e-52 1.444034e-53 1.440606e-54
## [6] 1.440264e-55 1.440230e-56 1.440226e-57 1.440226e-58 1.440226e-59
summary(abs(1 - y/1e-100))
plot(a, abs(1 - y/1e-100), log="xy", type="b")
stopifnot(abs(1 - y/1e-100) < 2e-13)# max (32b Linux, P.M) = 4.186e-14
## > The correct values of z1 should be:
z1true <- c(226.97154111939946, 222.15218724493326, 219.59373661415756,
217.27485383840451, 214.98015408183574, 212.68797118872064,
210.39614286838227, 208.10445550564617, 205.81289009100664,
203.52144711679352)
all.equal(z1[1:10], z1true, tol=0)# 1.307e-16 on 32-bit (Pentium M); 1.686e-16 (64b Lnx, 2019)
stopifnot(all.equal(z1[1:10], z1true, tol = 1e-15))
showProc.time()
##>
##> With these values of z1true, we get the expected values:
(y <- pgamma(z1true, a, lower.tail=FALSE))
## [1] 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100
## > I am using the precompiled binary version of R, under Windows Vista.
## > -----------
## >> version
## > _
## > platform i386-pc-mingw32
## > arch i386
## > os mingw32
## > system i386, mingw32
## > status
## > major 2
## > minor 7.1
## > year 2008
## > month 06
## > day 23
## > svn rev 45970
## > language R
## > version.string R version 2.7.1 (2008-06-23)
## > ------------
## >
## > So, it seems qgamma is inaccurate for small probability values in the upper
## > tail, when the shape parameter is also small.
###_-- MM: Still wrong:
(xm <- 2^-1074.9999) # is less than .Machine $ double.xmin == the really x > 0
pgamma(xm, .00001)# 0.992589
qgamma(.99, .00001)##--> NaN -- should give 0 or "xmin" or so
## FIXME -- ok, now
## but
curve(qgamma(x, .001, lower=FALSE), .4, .8, n=1001, log="y")
curve(qgamma(x, 1e-5, lower=FALSE), .002, .2, n=1001, log="xy")
curve(qgamma(x, 1e-7, lower=FALSE), 1e-5, .04, n=1001, log="xy")
curve(qgamma(x, 1e-12, lower=FALSE), 1e-12, 1e-2, n=1001, log="xy")
## or
curve(qgamma(x, 1e-121, lower=FALSE), 7e-119, 8e-119,
n=2001, log="y", yaxt="n")
try( # reveals eaxis() bug ? -- for the *subnormal* numbers
eaxis(2, at = 10^-seq(304,324, by=2))
)
curve(qgamma(x, .001, lower=FALSE), .4, .6, n=1001, log="y")
curve(qgamma(x, .001, lower=FALSE), .5, .55, n=1001, log="y")
## gives a *warning* from axis() because of subnormal y-range {was error, fixed in R-devel (2018-08)}
curve(qgamma(x, .001, lower=FALSE), .52, .53, n=1001, log="y")
curve(qgamma(x, .001, lower=FALSE)*1e100, .522, .526, n=1001, log="y")
showProc.time()
Try the DPQ package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.