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

In DPQ: Density, Probability, Quantile ('DPQ') Computations

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)

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.

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

Author(s)

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. \Sexpr[results=rd]{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")
})

DPQ documentation built on July 18, 2025, 3:01 p.m.