pbetaAS_eq20: `pbeta()` Approximations of Abramowitz & Stegun, $26.5.\{20,21\}$

Description

Compute (simple) pbeta() approximations; Abramowitz & Stegun, eq. 26.5.20 and 26.5.21. The first, eq. 26.5.20, approximates via a $\chi^2$ distribution, i.e., calls pchisq() with appropriate arguments, where the 2nd, eq. 26.5.21., approximates via the normal (aka “Gaussian”) distribution, i.e., calls pnorm().

In addition, rather for “didactical” reasons or simple comparisons, we also provide pbetaNorm2() which simply uses the normal approximation with matching moments, i.e., uses $Phi^{-1}((x - \mu)/\sigma)$ where $\mu = \mu(\alpha,\beta)$ is the mean, and $\sigma = \sigma(\alpha,\beta)$ is the standard deviation of the $B(\alpha,\beta)$ Beta distribution, where $\alpha=$shape1, $\beta=$shape2, $$\mu(\alpha,\beta) = \frac{\alpha}{\alpha+\beta},\ \ \mathrm{and}$$ $$\sigma(\alpha,\beta) = \frac{\sqrt{\alpha\beta},(\alpha+\beta)\sqrt{\alpha+\beta+1}}.$$ If the “validity” conditions are not fulfilled, the functions will try to “swap tails”, i.e., using the (anti-) symmetries of the Beta distribution, viz. computing $1 - F(1-x, b, a)$ instead of $F(x, a, b)$.

Usage

pbetaAS_eq20(q, shape1, shape2, lower.tail = TRUE, log.p = FALSE, warnBad = TRUE)
pbetaAS_eq21(q, shape1, shape2, lower.tail = TRUE, log.p = FALSE, warnBad = TRUE)
pbetaNorm2  (q, shape1, shape2, lower.tail = TRUE, log.p = FALSE)
## lower-level -- direct formula without range checks :
.pbeta.eq20(q, a,b, lower.tail = TRUE, log.p = FALSE, ._1_q = .5 - q + .5)
.pbetaSeq20(q, a,b, lower.tail = TRUE, log.p = FALSE)
.pbeta.eq21(q, a,b, lower.tail = TRUE, log.p = FALSE, ._1_q = .5 - q + .5)
.pbetaSeq21(q, a,b, lower.tail = TRUE, log.p = FALSE)

Value

a numeric vector, the length of e.g., length(q + shape1 + shape2), i.e., recycled to common length. With the corresponding pbeta() approximations, i.e., probabilities (in $[0,1]$) or log-probabilities in $[-\infty, 0]$.

Arguments

q, shape1, shape2: non-negative vectors (recycled to common length), q in $[0,1]$, see pbeta.
a, b: the same as shape1 and shape2, respectively; for convenience in (R) formulas only.
lower.tail: indicating if $F(q; *)$ should be returned or the upper tail probability $1 - F(q)$.
log.p: logical indicating if log(<pr>), log-probabilities should be returned instead of (default) the probabilities <pr>.
warnBad: logical indicating if a warning should be signalled when for some (q, shape1, shape2, lower.tail) combinations, A.&S. indicate the approximation to be bad, i.e., potentially quite inaccurate.
._1_q: $ = 1 - q$, computed numerically carefully (potentially with less cancellation error), or pre-computed.

Author

Martin Maechler

References

Abramowitz, M. and Stegun, I. A. (1972) Handbook of Mathematical Functions. New York: Dover. https://en.wikipedia.org/wiki/Abramowitz_and_Stegun provides links to the full text which is in public domain; formulas (26.5.20) & (26.5.21), p.945.

Beta Distribution (in Wikipedia) https://en.wikipedia.org/wiki/Beta_distribution

Examples

Run this code

str(alf <- 2^seq(1,15, by=1/8)) # 2 1.18 ... 32768
bet <- seq_along(alf) # 1 2 .. 113

E <- alf/(alf + bet) # = E[X]
pb   <- pbeta       (E, alf, bet)
pb20 <- pbetaAS_eq20(E, alf, bet) # warning .. 112 potentially bad
pb21 <- pbetaAS_eq21(E, alf, bet)
pbN2 <- pbetaNorm2  (E, alf, bet); summary(pbN2)# all == 0.5 as q = E[ Beta(a,b) ]
pbM <- cbind(pb20, pb21, pbN2)
tit <- quote(list(pbeta[...](E, alpha, beta) ~~" ", alpha == 2^{1~".."~15},
                  beta==1:113, E == frac(alpha,alpha+beta)))
matplot(alf, cbind(pb, pbM), type = "l", lwd = 1:4, log = "x",
        xlab = quote(alpha ~~ "[log]"), xaxt = "n", main = tit)
sfsmisc::eaxis(1, sub10=2)
legend("topright", ncol=2, bty="n", lty = 1:4, col = adjustcolor(1:4, 3/4), lwd = 1:4,
       expression(pbeta(), pbetaAS_eq20(), pbetaAS_eq21(), pbetaNorm2()))
stopifnot(apply(pbM, 2, all.equal, target = pb, tolerance = 0.03))
## relative errors, compared to R's pbeta()
round(cbind(sort(apply(pbM, 2, sfsmisc::relErr, target = pb))), 5)
## pb21 0.00013 << good here
## pb20 0.00356 << "ok" (got warning about "bad")
## pbN2 0.02742 << (constant)

## 2.  below E[]: mu - 2 s =======================================================
## -----------------------  --> use  xmin :
sdB <- sqrt(alf * bet)/(alf + bet)/sqrt(alf + bet + 1)
xmin <- unique(pmax(0, E - 2*sdB))
xmax <- unique(pmin(1, E + 2*sdB))
##
pb   <- pbeta       (xmin, alf, bet)
pb20 <- pbetaAS_eq20(xmin, alf, bet) # warning .. 104 potentially bad
## --> analysis shows the 'swap' logical to be clearly sub optimal.  ==> functions .pbeta.eq20() etc
pb21 <- pbetaAS_eq21(xmin, alf, bet)
pbN2 <- pbetaNorm2  (xmin, alf, bet); summary(pbN2)
pbM <- cbind(pb20, pb21, pbN2)
round(cbind(sort(apply(pbM, 2, sfsmisc::relErr, target = pb))), 5)
## pb21 0.00594
## pbN2 0.27069
## pb20 0.68668

tit <- quote(list(pbeta[...](xmin, alpha, beta) ~~" ", alpha == 2^{1~".."~15},
                  beta==1:113, xmin == (mu - 2*sigma)(alpha,beta)))
matplot(alf, cbind(pb, pbM), type = "l", lwd = 1:4, log = "x",
        xlab = quote(alpha ~~ "[log]"), xaxt = "n", main = tit)
sfsmisc::eaxis(1, sub10=2)
legend("topright", ncol=2, bty="n", lty = 1:4, lwd = 1:4, col = 1:4, # col = adjustcolor(1:4, 3/4),
       legend = expression(pbeta(), pbetaAS_eq20(), pbetaAS_eq21(), pbetaNorm2()))
##
pb20. <- .pbeta.eq20(xmin, alf,bet)
pb20S <- .pbetaSeq20(xmin, alf,bet)
pb21. <- .pbeta.eq21(xmin, alf,bet)
pb21S <- .pbetaSeq21(xmin, alf,bet)
matlines(alf, cbind(pb20., pb20S, pb21., pb21S), lwd = 4, lty = 1, col = adjustcolor(2:5, 1/4))
legend("right", ncol=2, bty="n", lty = 1, lwd = 4, col = adjustcolor(2:5, 1/4),
       legend = paste0(".pbeta", c(".eq20", "Seq20", ".eq21", "Seq21"), "()"))

## 3.  above E[]: mu + 2 s =======================================================
## -----------------------  --> use  xmax
##
pb   <- pbeta       (xmax, alf, bet)
pb20 <- pbetaAS_eq20(xmax, alf, bet) # warning .. 104 potentially bad
## --> analysis shows the 'swap' logical to be clearly sub optimal.  ==> functions .pbeta.eq20() etc
pb21 <- pbetaAS_eq21(xmax, alf, bet)
pbN2 <- pbetaNorm2  (xmax, alf, bet); summary(pbN2)
pbM <- cbind(pb20, pb21, pbN2)
round(cbind(sort(apply(pbM, 2, sfsmisc::relErr, target = pb))), 5)
## pb21 0.00009
## pb20 0.00559
## pbN2 0.00610
tit <- quote(list(pbeta[...](xmax, alpha, beta) ~~" ", alpha == 2^{1~".."~15},
                  beta==1:113, xmax == (mu + 2*sigma)(alpha,beta)))
matplot(alf, cbind(pb, pbM), type = "l", lwd = 1:4, log = "x",
        xlab = quote(alpha ~~ "[log]"), xaxt = "n", main = tit)
sfsmisc::eaxis(1, sub10=2)
legend("topright", ncol=2, bty="n", lty = 1:4, lwd = 1:4, col = 1:4, # col = adjustcolor(1:4, 3/4),
       legend = expression(pbeta(), pbetaAS_eq20(), pbetaAS_eq21(), pbetaNorm2()))
pb20. <- .pbeta.eq20(xmax, alf,bet)
pb20S <- .pbetaSeq20(xmax, alf,bet)
pb21. <- .pbeta.eq21(xmax, alf,bet)
pb21S <- .pbetaSeq21(xmax, alf,bet)
matlines(alf, cbind(pb20., pb20S, pb21., pb21S), lwd = 4, lty = 1, col = adjustcolor(2:5, 1/4))
legend("right", ncol=2, bty="n", lty = 1, lwd = 4, col = adjustcolor(2:5, 1/4),
       legend = paste0(".pbeta", c(".eq20", "Seq20", ".eq21", "Seq21"), "()"))

Run the code above in your browser using DataLab

pbetaAS_eq20: `pbeta()` Approximations of Abramowitz & Stegun, \(26.5.\{20,21\}\)