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

Description

Computes $$algdiv(a,b) := \log \frac{\Gamma(b)}{\Gamma(a+b)} = \log \Gamma(b) - \log\Gamma(a+b) = \code{lgamma(b) - lgamma(a+b)}$$ in a numerically stable way.

This is an auxiliary function ins TOMS 708 implementation of pbeta, aka the incomplete beta function ratio.

Usage

algdiv(a, b)

Arguments

a, b

numeric vectors which will be recycled to the same length.

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

Details

Note that this is also useful to compute the Beta function $$B(a,b) = \frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}.$$ Clearly, $$\log B(a,b) = \log\Gamma(a) + algdiv(a,b) = \log\Gamma(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() ?

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.

Examples

Run this code

# NOT RUN {
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)
# }

Run the code above in your browser using DataCamp Workspace