DPQ (version 0.5-9)

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

Description

Computes $$\code{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 in R's (TOMS 708) implementation of pbeta(), aka the incomplete beta function ratio.

Usage

algdiv(a, b)

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

Arguments

a, b

numeric vectors which will be recycled to the same length.

Author

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

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) + \mathrm{algdiv(a,b)} = \log\Gamma(a) - \mathrm{logQab}(a,b)$$

In our ../tests/qbeta-dist.R we look into computing \(\log(p B(p,q))\) accurately for \(p \ll q\) .

We are proposing a nice solution there.
How is this related to algdiv() ?

Additionally, we have defined $$Qab = Q_{a,b} := \frac{\Gamma(a+b),\Gamma(b)},$$ such that \(\code{logQab(a,b)} := \log Qab(a,b)\) fulfills simply $$\code{logQab(a,b)} = - \code{algdiv(a,b)}$$ see logQab_asy.

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), tolerance=0) # 1.283e-16 !!
(lQab <- logQab_asy(3, 1e10))
## relative error
1 + lQab/ algdiv(3, 1e10) # 0 (64b F 30 Linux; 2019-08-15)

## in-and outside of "certified" argument range {b >= 8}:
a. <- c(1:3, 4*(1:8))/32
b. <- seq(1/4, 20, by=1/4)
ad <- t(outer(a., b., algdiv))
## direct computation:
f.algdiv <- function(a,b) lgamma(b) - lgamma(a+b)
ad.d <- t(outer(a., b., f.algdiv))

matplot (b., ad.d, type = "o", cex=3/4,
         main = quote(log(Gamma(b)/Gamma(a+b)) ~"  vs.  algdiv(a,b)"))
mtext(paste0("a[1:",length(a.),"] = ",
        paste0(paste(head(paste0(formatC(a.*32), "/32")), collapse=", "), ", .., 1")))
matlines(b., ad,   type = "l", lwd=4, lty=1, col=adjustcolor(1:6, 1/2))
abline(v=1, lty=3, col="midnightblue")
# The larger 'b', the more accurate the direct formula wrt algdiv()
all.equal(ad[b. >= 1,], ad.d[b. >= 1,]       )# 1.5e-5
all.equal(ad[b. >= 2,], ad.d[b. >= 2,], tol=0)# 3.9e-9
all.equal(ad[b. >= 4,], ad.d[b. >= 4,], tol=0)# 4.6e-13
all.equal(ad[b. >= 6,], ad.d[b. >= 6,], tol=0)# 3.0e-15
all.equal(ad[b. >= 8,], ad.d[b. >= 8,], tol=0)# 2.5e-15 (not much better)