algdiv {DPQ}R Documentation

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

a, b

numeric vectors which will be recycled to the same length.

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)

[Package DPQ version 0.4-4 Index]