pnormLU {DPQ} | 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 vector. |
lower.tail, log.p |
logical, see, e.g., |
a numeric vector like x
Martin Maechler
Lutz Duembgen (2010)
Bounding Standard Gaussian Tail Probabilities;
arXiv preprint 1012.2063
,
https://arxiv.org/abs/1012.2063
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) 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)) col4 <- adjustcolor(1:4, 1/2) doLegTit <- function() { 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=col4, 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)) ## 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: x <- - 2^seq(0, 30, by = 1/64) ## == 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))