pp <- c(.001, .005, .01, .05, (1:9)/10, .95, .99, .995, .999)
z_p <- qnorm(pp)
assertDeprecation <- function(expr, verbose=TRUE)
tools::assertCondition(expr, verbose=verbose, "deprecatedWarning")
assertDeprecation(qA <- qnormAppr(pp))
(R <- cbind(pp, z_p, qA,
qUA = qnormUappr(pp, lower.tail= TRUE),
qA6 = qnormUappr6(pp, lower.tail=TRUE)))
## Errors, absolute and relative:
relEr <- function(targ, curr) { ## simplistic "smart" rel.error
E <- curr - targ
r <- E/targ # simple, but fix 0/0:
r[targ == 0 & E == 0] <- 0
r
}
mER <- cbind(pp,
errA = z_p - R[,"qA" ],
errUA = z_p - R[,"qUA"],
rE.A = relEr(z_p, R[,"qA" ]),
rE.UA = relEr(z_p, R[,"qUA"]),
rE.A6 = relEr(z_p, R[,"qA6"]))
signif(mER)
lp <- -c(1000, 500, 200, 100, 50, 20:10, seq(9.75, 0, by = -1/8))
signif(digits=5, cbind(lp # 'p' need not be specified if 'lp' is !
, p. = -expm1(lp)
, qnU = qnormUappr (lp=lp)
, qnU6= qnormUappr6(lp=lp)
, qnA1= qnormAsymp(lp=lp, lower.tail=FALSE, order=1)
, qnA5= qnormAsymp(lp=lp, lower.tail=FALSE, order=5)
, qn = qnorm(lp, log.p=TRUE)
)) ## oops! shows *BUG* for last values where qnorm() > 0 !
curve(qnorm(x, lower.tail=FALSE), n=1001)
curve(qnormUappr(x), add=TRUE, n=1001, col = adjustcolor("red", 1/2))
## Error curve:
curve(qnormUappr(x) - qnorm(x, lower.tail=FALSE), n=1001,
main = "Absolute Error of qnormUappr(x)")
abline(h=0, v=1/2, lty=2, col="gray")
curve(qnormUappr(x) / qnorm(x, lower.tail=FALSE) - 1, n=1001,
main = "Relative Error of qnormUappr(x)")
abline(h=0, v=1/2, lty=2, col="gray")
curve(qnormUappr(lp=x) / qnorm(x, log.p=TRUE) - 1, -200, -1, n=1001,
main = "Relative Error of qnormUappr(lp=x)"); mtext(" & qnormUappr6() [log.p scale]", col=2)
curve(qnormUappr6(lp=x) / qnorm(x, log.p=TRUE) - 1, add=TRUE, col=2, n=1001)
abline(h=0, lty=2, col="gray")
curve(qnormUappr(lp=x) / qnorm(x, log.p=TRUE) - 1,
-2000, -.1, ylim = c(-2e-4, 1e-4), n=1001,
main = "Relative Error of qnormUappr(lp=x)"); mtext(" & qnormUappr6() [log.p scale]", col=2)
curve(qnormUappr6(lp=x) / qnorm(x, log.p=TRUE) - 1, add=TRUE, col=2, n=1001)
abline(h=0, lty=2, col="gray")
## zoom out much more - switch x-axis {use '-x'} and log-scale:
curve(qnormUappr6(lp=-x) / qnorm(-x, log.p=TRUE) - 1,
.1, 1.1e10, log = "x", ylim = 2.2e-4*c(-2,1), n=2048,
main = "Relative Error of qnormUappr6(lp = -x) [log.p scale]") -> xy.q
abline(h=0, lty=2, col="gray")
## 2023-02: qnormUappr6() can be complemented with
## an approximation around center p=1/2: qnormCappr()
p <- seq(0,1, by=2^-10)
M <- cbind(p, qn=(qn <- qnorm(p)),
reC1 = relEr(qn, qnormCappr(p)),
reC2 = relEr(qn, qnormCappr(p, k=2)),
reC3 = relEr(qn, qnormCappr(p, k=3)),
reU6 = relEr(qn, qnormUappr6(p,lower.tail=TRUE)))
matplot(M[,"p"], M[,-(1:2)], type="l", col=2:7, lty=1, lwd=2,
ylim = c(-.004, +1e-4), xlab=quote(p), ylab = "relErr")
abline(h=0, col="gray", lty=2)
oo <- options(width=99)
summary( M[,-(1:2)])
summary(abs(M[,-(1:2)]))
options(oo)
Run the code above in your browser using DataLab