[R] The gamma function and infinity

[Martin Maechler]

>>>>> "PD" == Peter Dalgaard <p.dalgaard at biostat.ku.dk>
>>>>>     on 05 Dec 2005 17:14:47 +0100 writes:

    PD> Florent Bresson <f_bresson at yahoo.fr> writes:
    >> I have to calculate some formula like:
    >> 
    >> gamma(x)/(gamma(x+y)
    >> 
    >> and I observed that for relatively big values of x, R
    >> returns infinity and so cannot compute the formula. Is it
    >> possible to force R to give the real value of gamma(x)
    >> instead of Inf ?

    PD> No, it is overflowing the floating point representation.

    PD> How about using lgamma? That's basically what it is for.

Indeed.  Note however that in the case of really large x and
moderate y, even in the lgamma() formulation, there will be loss
of precision due to cancellation. There are "famous" asymptotic
(non-convergent series) formulae, you can use in that case, the
most simple one being
    Gamma(a + y) / Gamma(a) ~ a^y    (for a -> Inf,  y << a)
and hence
    log(Gamma(a + y)) - log(Gamma(a)) ~  y*ln(a)

--
Martin Maechler, ETH Zurich