DPQ (version 0.6-0)

dltgammaInc: TOMS 1006 - Fast and Accurate Generalized Incomplete Gamma Function

Description

This is the deltagammainc algorithm as described in Abergel and Moisan (2020).

It uses one of three different algorithms to compute \(G(p,x)\), their conveniently scaled generalization of the (upper and lower) incomplete Gamma functions,

$$G(p,x) = e^{x - p\log x} \times \gamma(p,x) \ \textrm{if} x \le p or \Gamma(p,x) \textrm{otherwise},$$

for \(\gamma(p,x)\) and \(\Gamma(p,x)\) the lower and upper incomplete gamma functions, and for all \(x \ge 0\) and \(p > 0\).

dltgammaInc(x,y, mu, p) := \(I_{x,y}^{\mu,p}\) computes their (generalized incomplete gamma) integral $$I_{x,y}^{\mu,p} := \int_x^y s^{p-1} e^{-\mu s} \; ds.$$

Current explorations (via package Rmpfr's igamma()) show that this algorithm is less accurate than R's own pgamma() at least for small \(p\).

Usage

dltgammaInc(x, y, mu = 1, p)

Value

a list with three components

rho, sigma

numeric vectors of same length (the common length of x,y,.., after recycling). ans := rho * exp(sigma) are the computed values of \(G(p,x)\).

method

integer code indicating the algorithm used for the computation of (rho, sigma).

Arguments

x,y, p

numeric vectors, typically of the same length; otherwise recycled to common length.

mu

a number different from 0; must be finite; 1 by default.

Author

C source from Remy Abergel and Lionel Moisan, see the reference; original is in 1006.zip at https://calgo.acm.org/.

Interface to R in DPQ: Martin Maechler

References

Rémy Abergel and Lionel Moisan (2020) Algorithm 1006: Fast and accurate evaluation of a generalized incomplete gamma function, ACM Transactions on Mathematical Software 46(1): 1--24. tools:::Rd_expr_doi("10.1145/3365983")

See Also

Rmpfr's igamma(), R's pgamma().

Examples

p <- pi
x <- y <- seq(0, 10, by = 1/8)

dltg1 <- dltgammaInc(0, y,   mu = 1, p = pi); r1 <- with(dltg1, rho * exp(sigma))  ##
dltg2 <- dltgammaInc(x, Inf, mu = 1, p = pi); r2 <- with(dltg2, rho * exp(sigma))
str(dltg1)
stopifnot(identical(c(rho = "double", sigma = "double", method = "integer"),
                    sapply(dltg1, typeof)))
summary(as.data.frame(dltg1))
pg1 <- pgamma(q = y, shape = pi, lower.tail=TRUE)
pg2 <- pgamma(q = x, shape = pi, lower.tail=FALSE)

stopifnot(all.equal(r1 / gamma(pi), pg1, tolerance = 1e-14)) # seen 5.26e-15
stopifnot(all.equal(r2 / gamma(pi), pg2, tolerance = 1e-14)) #  "   5.48e-15

if(requireNamespace("Rmpfr")) withAutoprint({
    asN  <- Rmpfr::asNumeric
    mpfr <- Rmpfr::mpfr
    iGam <- Rmpfr::igamma(a= mpfr(pi, 128),
                          x= mpfr(y,  128))
    relErrV <- sfsmisc::relErrV
    relE <- asN(relErrV(target = iGam, current = r2))
    summary(relE)
    plot(x, relE, type = "b", col=2, # not very good: at  x ~= 3, goes up to  1e-14 !
         main = "rel.Error(incGamma(pi, x)) vs  x")
    abline(h = 0, col = adjustcolor("gray", 0.75), lty = 2)
    lines(x, asN(relErrV(iGam, pg2 * gamma(pi))), col = adjustcolor(4, 1/2), lwd=2)
    ## ## ==> R's pgamma() is *much* more accurate !!  ==> how embarrassing for TOMS 1006 !?!?
    legend("topright", expression(dltgammaInc(x, Inf, mu == 1, p == pi),
                                  pgamma(x, pi, lower.tail== F) %*% Gamma(pi)), 
           col = c(palette()[2], adjustcolor(4, 1/2)), lwd = c(1, 2), pch = c(1, NA), bty = "n")
})