Compute log(gamma(b)/gamma(a+b)) when b >= 8

Description

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.

Usage

algdiv(a, b)

Arguments

Details

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() ?
 

Value

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).

Author(s)

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

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.

See Also

gamma, beta; my own logQab_asy().

Examples

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)