Nothing
#### Extreme-tail pnorm() and qnorm() approximations
#### ===============================================
### notably pnormAsymp() and qnormAsymp()
### ~~~~~~~~~~~~ ~~~~~~~~~~~~
## partly inspired by MM's NUMERICS/dpq-functions/qnorm-asymptotic.R
## and also ../man/pnormAsymp.Rd
## ../man/qnormAsymp.Rd
library(DPQ)
library(sfsmisc)# {we import it in DPQ}
p.ver <- function() mtext(R.version.string, cex=3/4, adj=1)
readUser <- function(i, max.i) {
if(interactive() && i != max.i) {
cat("[Enter] to continue: ")
cat(readLines(stdin(), n=1), "\n")
}
}
# "bisque4" , "coral2" ..
p.epsC <- function(col = "coral2", alpha = 3/4, lwd1 = 3) {
lg <- par("ylog")
mult <- if(lg) c(1/2, 1, 2) else -2:2 # multiplicities
is1 <- abs(mult) == 1
abline(h = .Machine$double.eps * mult,
col = adjustcolor(col, alpha),
lty = ifelse(is1, 2, 5),
lwd = ifelse(is1, lwd1, 1))
}
absP <- function(re) pmax(abs(re), 2e-17) # not zero, so log-scale "shows" it
if(!dev.interactive(orNone=TRUE)) pdf("pqnorm-extreme.pdf")
x <- sort(2 ^ c(seq(3, 10, by=1/16), seq(10.125, 17, by=1/8)))
str(x)
summary(x)
pn <- pnorm(x, lower.tail=FALSE, log.p=TRUE)
## look at the asymptotic approximations:
k.s <- 0:5; nks <- paste0("k=", k.s); k. <- setNames(k.s, nks)
pnAsy <- sapply(k., function(k) pnormAsymp(x, k=k, lower.tail=FALSE, log.p=TRUE))
relEpnAsy <- apply(pnAsy, 2, relErrV, target = pn)
tit <- "relErrV(pnormAsym(x, k=k, lower.tail=FALSE, log.p=TRUE)\n vs pnorm(x, ..)"
leg.k <- function(x, y=NULL, ...) legend(x,y, legend=paste0("k=",0:5), lwd=2, col=1:6, lty=1:6, bty="n", ...)
matpl <- function(x, y, ...) {
matplot(x,y, type="l", lwd=2, log="x", main=tit, xaxt="n", ...)
eaxis(1, sub10=3); grid()
}
matpl(x, relEpnAsy, ylim = c(-8,1)*5e-6 ); leg.k("center")
matpl(x, relEpnAsy, ylim = c(-1,1)*1e-8 ); leg.k("topright")
matpl(x, relEpnAsy, ylim = c(-1,1)*1e-10); leg.k("topright")
matpl(x, relEpnAsy, ylim = c(-1,1)*1e-12); leg.k("topright")
matpl(x, relEpnAsy, ylim = c(-1,1)*1e-14); leg.k("topright"); p.epsC()
matpl(x, relEpnAsy, ylim = c(-1,1)*1e-15); leg.k("topright"); p.epsC()
abline(v = 38581.371, col=adjustcolor(4, 1/2), lwd=3) # the Rmpfr maximal value (for "independent accurate pnorm")
## They seem to be of opposite sign - mostly (and apart from +/- epsC)
## *but* k=5 (at the beginning; (and then +/- epsC)
## ?proof? --- use Rmpfr-version pnormAsymp() but there,
## ....... pnorm(x, log.p=TRUE) *underflows* to -Inf as soon as x > 38581.371
## Start with quantiles so we know the "truth according to pnorm()" :
qs <- c(2^seq(0, 35, by=1/256), Inf) # => s >= 1.84 --> no NaN in xs_5 etc
## == the "true" qnorm() according to pnorm()
lp <- pnorm(qs, lower.tail=FALSE, log.p=TRUE)
qnp <- qnorm(lp, lower.tail=FALSE, log.p=TRUE)
s <- -lp # = -log(1 - Phi(qs))
summary(r <- sqrt(s))
## store in 'dat' including current "inaccuracy" of qnorm():
head(dat <- data.frame(l2q=log2(qs), qs, # lp,
s, # r = sqrt(s), t = log(s),
relE_qn = relErrV(qs, qnp)))
##^^^^^^ will depend much on R version
stopifnot(dat[nrow(dat), "relE_qn"] == 0)
## relative Error (log x)
plot(relE_qn ~ s, dat, log="x", type="l", col = 2); p.ver()
## |rel.Error| ==> log-log scale
plot(abs(relE_qn) ~ s, dat, log="xy", type="l", col = 2,
axes=FALSE, main = "|relative Error| of qnorm(-s, log.p=TRUE, lower=F) wrt pnorm()")
eaxis(1); eaxis(2); p.ver(); p.epsC()
## now look at the asymptotic approximations:
k.s <- 0:5; nks <- paste0("k=", k.s); k. <- setNames(k.s, nks)
qnAsy <- sapply(k., function(ord) qnormAsymp(lp=lp, order=ord))
stopifnot(identical(qnAsy,
sapply(k., function(ord) qnormAsymp(p=lp, lower.tail=FALSE, log.p=TRUE, order=ord))
))
relEAsy <- apply(qnAsy, 2, relErrV, target = qs)
matplot(sqrt(s), relEAsy, type="l", log="x", xlab = quote(r == sqrt(s)),
main = "relative Error of qnormAsymp(.., k=*)"); p.ver()
legend("top", nks, col=1:6, lty=1:5, bty="n")
mp <- c(-1,1)
xL <- list(NULL, c(2, 1000), c(2, 1e5), c(5, 1e5), c(10, 1e8), c(10, 1e10))
yL <- list(NULL, mp * 0.1, mp * 1e-4, mp * 1e-7, mp * 1e-11, mp * 2e-15)
stopifnot(length(xL) == length(yL))
for(j in seq_along(xL)) {
matplot(sqrt(s), relEAsy, type="l", log="x", xaxt="n", lwd=2,
xlab = quote(r == sqrt(s)),
xlim = xL[[j]], ylim = yL[[j]],
main = "relative Error of qnormAsymp(.., k=*)")
eaxis(1, sub10=2) # not here: p.ver()
legend("top", nks, col=1:6, lty=1:5, lwd=2, bty="n")
readUser(j, length(xL))
}
for(j in seq_along(xL)) {
if(!is.null(yli <- yL[[j]])) yli <- c(2e-17, 100*yli[2])
if(!is.null(xli <- xL[[j]])) xli[2] <- 10*xli[2]
cat("j=",j,"; rbind(xlim,ylim) = \n"); print(rbind(xlim=xli, ylim=yli))
cols <- adjustcolor(1:6, 3/4)
## abs() does not show 0 {-> -Inf); absP() will show them
matplot(sqrt(s), abs(relEAsy), type="l", log="xy", xaxt="n", yaxt="n", col=cols, lwd=2,
xlab = quote(r == sqrt(s)),
xlim = xli, ylim = yli,
main = "relative Error of qnormAsymp(.., k=*)")
matlines(sqrt(s), absP(relEAsy), col=adjustcolor(cols, 1/2))
eaxis(1, sub10=2) # not here: p.ver()
eaxis(2, sub10=c(-3,2))
legend("top", nks, col=cols, lwd=2, lty=1:5, bty="n")
readUser(j, length(xL))
}
## Real tests :
table(relEAsy[s > 1e17,] * 2^52)
## -0.5 0 1
## 71 9434 125
table(relEAsy[s > 1e16,] * 2^52)
## -0.5 0 1 2 3 4 5
## 71 11550 311 135 78 18 17
## drop k = 0 from now on ------------- The following *is* platform dependent .. let's see
re1 <- abs(relEAsy[s > 95e6, "k=1"])
table(re1 * 2^52)
## 0 0.5 1
## 5264 84 93
plot(abs(relEAsy[,"k=2"]) ~ s, subset = s > 1e5 & s < 1e8, log="xy", type="l", xaxt="n")
eaxis(1); p.epsC()
re2 <- abs(relEAsy[s > 2e5, "k=2"])
table(re2 * 2^52)
## 0 0.5 1
## 6353 85 141
plot(abs(relEAsy[,"k=3"]) ~ s, subset = s > 500 & s < 1e5, log="xy", type="l", xaxt="n")
eaxis(1); p.epsC()
re3 <- abs(relEAsy[s > 4000, "k=3"])
table(re3 * 2^52)
## 0 0.5 1
## 6953 124 225
plot(abs(relEAsy[,"k=4"]) ~ s, subset = s > 300 & s < 2000, log="xy", type="l", axes=FALSE)
eaxis(1);eaxis(2); p.epsC()
re4 <- abs(relEAsy[s > 1500, "k=4"])
table(re4 * 2^52)
## 0 0.5 1
## 7154 145 185
plot(abs(relEAsy[,"k=5"]) ~ s, subset = s > 200 & s < 1000, log="xy", type="l", axes=FALSE)
eaxis(1);eaxis(2); p.epsC()
r. <- pretty(sqrt( par("xaxp")[1:2] ), 10)
(rlab <- as.expression(lapply(r., function(rr) substitute(R^2, list(R = rr+0)))))
axis(3, at=r.^2, labels = rlab, col = 4, col.axis = 4, mgp = c(1.25,.5,0))
re5 <- abs(relEAsy[s > 700, "k=5"])
table(re5 * 2^52)
## 0 0.5 1
## 7285 141 199
stopifnot(exprs = {
abs(relEAsy[s > 1e17,]) <= 2^-52 # even for k=0
re1 <= 2^-52
re2 <= 2^-52
re3 <= 2^-52
re4 <= 2^-52
re5 <= 2^-52
})
### R code from vignette source '../vignettes/qnorm-asymp.Rnw' ====================
###
### code chunk number 2: qnormLog-compute
#########################################
qs <- 2^seq( 0, 29, by=1/256) # => s >= 1.84
lp <- pnorm(qs, lower.tail=FALSE, log.p=TRUE)
s <- -lp # = -pnorm(..) = -log(1 - Phi(qs)) > 0
##
qnrm <- qnorm (-s, lower.tail=FALSE, log.p=TRUE)
qnrm405 <- qnormR(-s, lower.tail=FALSE, log.p=TRUE, version= "4.0.x") # R <= 4.0.5
qnrm410 <- qnormR(-s, lower.tail=FALSE, log.p=TRUE, version= "2020-10-17")
Rver <- shortRversion()
if(getRversion() <= "4.0.5") { # our qnormR(.., version="4.0.x")
cat(sprintf("%s, \"4.0.5\",\n all.equal(*, tol=0): %s; identical(): %s\n", Rver,
all.equal(qnrm, qnrm405, tolerance=0), identical(qnrm, qnrm405)))
stopifnot(all.equal(qnrm, qnrm405, tolerance = 1e-12))
} else if(getRversion() < "4.3") { # our qnormR(*, version="2020-10-17") matches:
cat(sprintf("%s, \"4.1.0\",\n all.equal(*, tol=0): %s; identical(): %s\n", Rver,
all.equal(qnrm, qnrm410, tolerance=0), identical(qnrm, qnrm410)))
## see TRUE twice, for R 4.2.2, Linux{x86_64-pc-linux-gnu} *and* Windows{x86_64-w64-mingw32/x64}
# M1mac(aarch64-apple-darwin20, R 4.2.1 ptchd): 2.675587e-16
stopifnot(all.equal(qnrm, qnrm410, tolerance = 1e-12))
} else { # R version >= 4.3.x
qnrm43 <- qnormR(-s, lower.tail=FALSE, log.p=TRUE, version = "2022")
cat(sprintf("%s, >= 4.3.x,\n all.equal(*, tol=0): %s; identical(): %s\n", Rver,
all.equal(qnrm, qnrm43, tolerance=0), identical(qnrm, qnrm43)))
rE6 <- qnorm(-1e6, log.p=TRUE)/-1414.2077829910174 - 1
cat(sprintf(" rE(-1e6) = %g\n", rE6))
if(abs(rE6) < 7e-16) # have R-devel with new 2022 code:
stopifnot(all.equal(qnrm, qnrm43, tolerance = 1e-14))
}
source(system.file("extraR", "qnorm-asymp-utils.R", package="DPQ"))
if(!exists("r0", mode="numeric"))
q("no")
if(!dev.interactive(orNone=TRUE)) {
dev.off()
pdf("pqnorm-qnormAsy2.pdf", height=8*sqrt(2), width=8) # ~ 'A4'
}
###################################################
### code chunk number 27: plot-qnormAsy2
###################################################
sfsmisc::mult.fig(5, main = "qnormAsymp(*, k) approximations in the 5 cutpoint regions")
r0 <- c(27, 55, 109, 840, 36000, 6.4e8) # <-- cutoffs <--> in ../R/norm_f.R
# use k = 5 4 3 2 1 0 e.g. k = 0 good for r >= 6.4e8
for(ir in 2:length(r0))
p.qnormAsy2(r0[ir], k = 5 +2-ir, cex.main = .90)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.