pbeta_ser | R Documentation |
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.
pbeta_ser(q, shape1, shape2, log.p = FALSE, eps = 1e-15, errPb = 0, verbose = FALSE)
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. |
pbeta_ser()
crucially needs three auxiliary functions which we
“mpfr-ized” as well: gam1M()
,
lgamma1pM()
, and algdivM
.
An approximation to the Beta probability P[X \le q]
for X \sim B(a,b),
(where a=
shape1
, and b=
shape2
).
Didonato and Morris and R Core team; separate packaging by Martin Maechler.
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; \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1145/131766.131776")}.
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.
(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))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.