Re: Problems with gamma, lgamma, and alternatives

[Mike Miller]

On Tue, 25 Mar 2003, Przemek Klosowski wrote:

>    octave:1> gamma(1000)
>            ...
>     SLATEC     DGAMMA     X SO BIG GAMMA OVERF         3         2         1
>
>    These usually occur during intermediate iterations of a minimization 
> routine,
>    where it's difficult to control the value that will be passed to gamma. The
>
> It's a tough one: on one hand, gamma(171) is largest double precision
> representable integer argument (~7.2574e+306), gamma(1000) would be, oh,
> about 999!, or ~4.02e+2564. If you are venturing there during a search
> you are hosed---there is no numerical accuracy there at all.


I'm curious.  I understand the problem with gamma and large values, but
what was the problem with lgamma?  Using my old version of Octave, I see
that it works for values up to about 10^305...

lgamma(10^305.4) =  1.7639e+308

After that it does fail in a similar way to gamma.  Do you really need to
take a gamma function of a number greater than 10^305?

If so, remember Stirling's approximation!  It will work well with large
numbers:

n! ~ sqrt(2*pi*n)*n^n*exp(-n)

Where '~' means "approximately equals".  Of course gamma(n) = (n-1)!.

Maybe Stirling's formula can be incorporated into Michael Creel's
algorithms to avoid problems when gamma is too large.  It depends on how
the gamma is being used.

It is also possible to develop an function llgamma(n) = log(lgamma(n))
that uses Stirling's approximation.  For n > 10^305 use

lgamma(n) = log(gamma(n)) ~ n*(log(n)-1)

and

llgamma(n) = log(log(gamma(n))) ~ log(n)+log(log(n)-1)

But might not help as much as I was hoping because these functions will
both return 'Inf' when n is sufficiently large.  To use llgamma
efficiently with giant numbers, you'd have to first transform to log(n):

x = log(n)
f(x) = x + log(x + 1)

But log(n) will return 'Inf' so you have to parse the number somehow
before taking the log.  E.g., to take the llgamma of 6.384e+4356 we would
do this:

octave:1> x=log(6.384)+4356*log(10);
octave:2> x+log(x+1)
ans =  1.0041e+04

Thus gamma(6.384e+4356) = exp(exp(1.0041e+04)).  If it would be more
useful, this could all be done in log10 instead.

Also note that lgamma and llgamma are not useful for computing ratios of
gammas.  I'm not really sure what they would be used for, but I assume
there is a reason to compute gamma functions of monsterous numbers!

Best,

Mike



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------