algdiv {DPQ} | R Documentation |
Computes
algdiv(a,b) := \log \frac{Γ(b)}{Γ(a+b)} = \log Γ(b) - \logΓ(a+b) = \code{lgamma(b) - lgamma(a+b)}
in a numerically stable way.
This is an auxiliary function in R's (TOMS 708) implementation of
pbeta()
, aka the incomplete beta function ratio.
algdiv(a, b)
a, b |
numeric vectors which will be recycled to the same length. |
Note that this is also useful to compute the Beta function
B(a,b) = Γ(a)Γ(b)/Γ(a+b).
Clearly,
\log B(a,b) = \logΓ(a) + algdiv(a,b) = \logΓ(a) - logQab(a,b)
In our ../tests/qbeta-dist.R we look into computing log(p*Beta(p,q)) accurately for p << q --------------------- We are proposing a nice solution there. How is this related to algdiv() ?
a numeric vector of length max(length(a), length(b))
(if neither
is of length 0, in which case the result has length 0 as well).
Didonato, A. and Morris, A., Jr, (1992); algdiv()
's C
version from the R sources, authored by the R core team; C and R
interface: 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.
gamma
, beta
;
my own logQab_asy()
.
Qab <- algdiv(2:3, 8:14) cbind(a = 2:3, b = 8:14, Qab) # recycling with a warning ## algdiv() and my logQab_asy() give *very* similar results for largish b: all.equal( - algdiv(3, 100), logQab_asy(3, 100), tol=0) # 1.283e-16 !! (lQab <- logQab_asy(3, 1e10)) ## relative error 1 + lQab/ algdiv(3, 1e10) # 0 (64b F 30 Linux; 2019-08-15)