pbeta_ser {DPQmpfr} | R Documentation |
Beta Distribution Function – ‘BPSER’ Series Expansion from TOMS 708
Description
Compute a version of the Beta cumulative distribution function
(pbeta()
in R), namely using the series expansion, named
BPSER()
, from “TOMS 708”, i.e., Didonato and Morris (1992).
This “pure R” function exists for didactical or documentational reasons on one hand,
as R's own pbeta()
uses this expansion when appropriate and
other algorithms otherwise.
On the other hand, using high precision q
and MPFR arithmetic (via
package Rmpfr) may allow to get highly accurate pbeta()
values.
Usage
pbeta_ser(q, shape1, shape2, log.p = FALSE, eps = 1e-15, errPb = 0, verbose = FALSE)
Arguments
q , shape1 , shape2 |
quantiles and shape parameters of the Beta
distribution, |
log.p |
if TRUE, probabilities |
eps |
non-negative number; |
errPb |
an integer code, typically in |
verbose |
logical indicating if console output about intermediate results should be printed. |
Details
pbeta_ser()
crucially needs three auxiliary functions which we
“mpfr-ized” as well: gam1M()
,
lgamma1pM()
, and algdivM
.
Value
An approximation to the Beta probability P[X \le q]
for X \sim B(a,b),
(where a=
shape1
, and b=
shape2
).
Author(s)
Didonato and Morris and R Core team; separate packaging by Martin Maechler.
References
Didonato, A. and Morris, A., Jr, (1992) Algorithm 708: Significant digit computation of the incomplete beta function ratios, ACM Transactions on Mathematical Software 18, 360–373; doi:10.1145/131766.131776.
See Also
pbeta
, DPQmpfr's own pbetaD94
;
even more pbeta()
approximations in package DPQ, e.g.,
pnbetaAS310
, or pbetaRv1
.
In addition, for integer shape parameters, the potentially “fully accurate”
finite sum base pbetaI()
in package Rmpfr.
Examples
(p. <- pbeta_ser(1/2, shape1 = 2, shape2 = 3, verbose=TRUE))
(lp <- pbeta_ser(1/2, shape1 = 2, shape2 = 3, log.p = TRUE))
all.equal(lp, log(p.), tolerance=0) # 1.48e-16
stopifnot(all.equal(lp, log(p.), tolerance = 1e-13))
## Using Vectorize() in order to allow vector 'q' e.g. for curve():
str(pbetaSer <- Vectorize(pbeta_ser, "q"))
curve(pbetaSer(x, 1.5, 4.5)); abline(h=0:1, v=0:1, lty=2, col="gray")
curve(pbeta (x, 1.5, 4.5), add=TRUE, col = adjustcolor(2, 1/4), lwd=3)
## now using mpfr-numbers:
half <- 1/Rmpfr::mpfr(2, 256)
(p2 <- pbeta_ser(half, shape1 = 1, shape2 = 123))