pnormLU {DPQmpfr} | R Documentation |
Bounds for 1 - Φ(x), i.e., pnorm(x, *,
lower.tail=FALSE)
, typically related to Mill's Ratio.
pnormL_LD10(x, lower.tail = FALSE, log.p = FALSE) pnormU_S53 (x, lower.tail = FALSE, log.p = FALSE)
x |
positive (at least non-negative) numeric |
lower.tail, log.p |
logical, see, e.g., |
vector/array/mpfr like x
.
Martin Maechler
Lutz Duembgen (2010)
Bounding Standard Gaussian Tail Probabilities;
arXiv preprint 1012.2063
,
https://arxiv.org/abs/1012.2063
pnorm
. The same functions “numeric-only” are in my
DPQ package.
x <- seq(1/64, 10, by=1/64) px <- cbind( lQ = pnorm (x, lower.tail=FALSE, log.p=TRUE) , Lo = pnormL_LD10(x, lower.tail=FALSE, log.p=TRUE) , Up = pnormU_S53 (x, lower.tail=FALSE, log.p=TRUE)) matplot(x, px, type="l") # all on top of each other matplot(x, (D <- px[,2:3] - px[,1]), type="l") # the differences abline(h=0, lty=3, col=adjustcolor(1, 1/2)) ## check they are lower and upper bounds indeed : stopifnot(D[,"Lo"] < 0, D[,"Up"] > 0) matplot(x[x>4], D[x>4,], type="l") # the differences abline(h=0, lty=3, col=adjustcolor(1, 1/2)) ### zoom out to larger x : [1, 1000] x <- seq(1, 1000, by=1/4) px <- cbind( lQ = pnorm (x, lower.tail=FALSE, log.p=TRUE) , Lo = pnormL_LD10(x, lower.tail=FALSE, log.p=TRUE) , Up = pnormU_S53 (x, lower.tail=FALSE, log.p=TRUE)) matplot(x, px, type="l") # all on top of each other matplot(x, (D <- px[,2:3] - px[,1]), type="l") # the differences abline(h=0, lty=3, col=adjustcolor(1, 1/2)) ## check they are lower and upper bounds indeed : table(D[,"Lo"] < 0) # no longer always true table(D[,"Up"] > 0) ## not even when equality (where it's much better though): table(D[,"Lo"] <= 0) table(D[,"Up"] >= 0) ## *relative* differences: matplot(x, (rD <- 1 - px[,2:3] / px[,1]), type="l", log = "x") abline(h=0, lty=3, col=adjustcolor(1, 1/2)) ## abs() matplot(x, abs(rD), type="l", log = "xy", axes=FALSE, # NB: curves *cross* main = "relative differences 1 - pnormUL(x, *)/pnorm(x,*)") legend("top", c("Low.Bnd(D10)", "Upp.Bnd(S53)"), bty="n", col=1:2, lty=1:2) sfsmisc::eaxis(1, sub10 = 2) sfsmisc::eaxis(2) abline(h=(1:4)*2^-53, col=adjustcolor(1, 1/4)) ### zoom out to LARGE x : --------------------------- x <- 2^seq(0, 30, by = 1/64) col4 <- adjustcolor(1:4, 1/2) if(FALSE)## or even HUGE: x <- 2^seq(4, 513, by = 1/16) px <- cbind( lQ = pnorm (x, lower.tail=FALSE, log.p=TRUE) , a0 = dnorm(x, log=TRUE) , a1 = dnorm(x, log=TRUE) - log(x) , Lo = pnormL_LD10(x, lower.tail=FALSE, log.p=TRUE) , Up = pnormU_S53 (x, lower.tail=FALSE, log.p=TRUE)) doLegTit <- function(col=1:4) { title(main = "relative differences 1 - pnormUL(x, *)/pnorm(x,*)") legend("top", c("phi(x)", "phi(x)/x", "Low.Bnd(D10)", "Upp.Bnd(S53)"), bty="n", col=col, lty=1:4) } ## *relative* differences are relevant: matplot(x, (rD <- 1 - px[,-1] / px[,1]), type="l", log = "x", ylim = c(-1,1)/2^8, col=col4) ; doLegTit() abline(h=0, lty=3, col=adjustcolor(1, 1/2)) if(x[length(x)] > 1e150) # the "HUGE" case (not default) print( tail(cbind(x, px), 20) ) ##--> For very large x ~= 1e154, the approximations overflow *later* than pnorm() itself !! ## abs(rel.Diff) ---> can use log-log: matplot(x, abs(rD), type="l", log = "xy", xaxt="n", yaxt="n"); doLegTit() sfsmisc::eaxis(1, sub10=2) sfsmisc::eaxis(2) abline(h=(1:4)*2^-53, col=adjustcolor(1, 1/4)) ## lower.tail=TRUE (w/ log.p=TRUE) works "the same" for x < 0: require(Rmpfr) x <- - 2^seq(0, 30, by = 1/64) ## == log1mexp <- Rmpfr::log1mexp # Rmpfr version >= 0.8-2 (2020-11-11 on CRAN) px <- cbind( lQ = pnorm (x, lower.tail=TRUE, log.p=TRUE) , a0 = log1mexp(- dnorm(-x, log=TRUE)) , a1 = log1mexp(-(dnorm(-x, log=TRUE) - log(-x))) , Lo = log1mexp(-pnormL_LD10(-x, lower.tail=TRUE, log.p=TRUE)) , Up = log1mexp(-pnormU_S53 (-x, lower.tail=TRUE, log.p=TRUE)) ) matplot(-x, (rD <- 1 - px[,-1] / px[,1]), type="l", log = "x", ylim = c(-1,1)/2^8, col=col4) ; doLegTit() abline(h=0, lty=3, col=adjustcolor(1, 1/2)) ## Comparison with Rmpfr::erf() / erfc() based pnorm(): ## Set the exponential ranges to maximal -- to evade underflow as long as possible .mpfr_erange_set(value = (1-2^-52) * .mpfr_erange(c("min.emin","max.emax"))) l2t <- seq(0, 32, by=1/4) twos <- mpfr(2, 1024)^l2t Qt <- pnorm(twos, lower.tail=FALSE) pnU <- pnormU_S53 (twos, log.p=TRUE) pnL <- pnormL_LD10(twos, log.p=TRUE) logQt <- log(Qt) M <- cbind(twos, Qt, logQt = logQt, pnU) roundMpfr(M, 40) dM <- asNumeric(cbind(dU = pnU - logQt, dL = logQt - pnL, # NB: the numbers are *negative* rdU= 1 - pnU/logQt, rdL = pnL/logQt - 1)) data.frame(l2t, dM) ## The bounds are ok (where Qt does not underflow): L < p < U : stopifnot(pnU > pnL, pnU > logQt, (logQt > pnL)[Qt > 0]) roundMpfr(cbind(twos, pnL, pnU, D=pnU-pnL, relD=(pnU-pnL)/((pnU+pnL)/2)), 40) ## ----- R's pnorm() -- is it always inside [L, U] ?? --------------------- nQt <- stats::pnorm(asNumeric(twos), lower.tail=FALSE, log.p=TRUE) data.frame(l2t, check.names=FALSE , nQt , "L <= p" = c(" ", "W")[2 -(pnL <= nQt)] , "p <= U" = c(" ", "W")[2- (nQt <= pnU)]) ## ==> pnorm() is *outside* sometimes for l2t >= 7.25; always as soon as l2t >= 9.25 ## *but* the relative errors are around c_epsilon in all these cases : plot (2^l2t, asNumeric(abs(nQt-pnL)/abs(pnU)), type="o", cex=1/4, log="xy", axes=FALSE) sfsmisc::eaxis(1, sub10 = 2) sfsmisc::eaxis(2) lines(2^l2t, asNumeric(abs(nQt-pnU)/abs(pnU)), type="o", cex=1/4, col=2) abline(h=c(1:4)*2^-53, lty=2, col=adjustcolor(1, 1/4))