[SciPy-dev] Definition of gammaln(x) for negative x
G-J van Rooyen
gvrooyen@gmail....
Sat Nov 1 14:53:56 CDT 2008
I think it probably makes sense to keep it the way it is. AFAIK
gammaln is typically used to calculate products and divisions with the
gamma function, where it makes more sense to transform the entire
calculation into the log-domain, in order to prevent numerical
inaccuracies.
e.g. gamma(A)*gamma(B)/gamma(C)*gamma(D)
= exp(gammaln(A)+gammaln(B)-gammaln(C)+gammaln(D))
which works fine if all the arguments produce positive gamma-values.
If negative gamma-values are produced (as described in ticket #737),
the same calculation can be done by calculating gammaln on the
absolute value of the arguments, and doing a sign correction at the
end. For this, only the signs of gamma(A), etc. are needed. The
original scipy/special/cephes/gamma.c writes the sign to a global
variable named sgngam; this never gets imported to python.
It might make sense to keep gammaln as it is, but to optionally return
the sign information of gamma(A) in some way.
Your thoughts?
G-J
2008/11/1 David Cournapeau <david@ar.media.kyoto-u.ac.jp>:
> G-J van Rooyen wrote:
>> Hey everyone
>>
>> Ticket #737 refers:
>>
>> -----
>>
>> Gamma in negative for negative x whose floor value is odd. As such,
>> gammaln does not make sense for those values (while staying in the
>> real domain, at least). scipy.special.gammaln returns bogus values:
>>
>> import numpy as np
>> from scipy.special import gamma, gammaln
>> print np.log(gamma(-0.5))
>> print gammaln(-0.5)
>>
>> Returns nan in the first case (expected) and 1.26551212348 in the
>> second (totally meaningless value).
>>
>> -----
>>
>> The info line for gammaln reads:
>> * gammaln -- Log of the absolute value of the gamma function.
>>
>> With this definition of gammaln, the function actually works fine,
>> since np.log(abs(gamma(-0.5))) is in fact 1.2655.
>
> I have just checked with R, R does define log gamma as the
> log(abs(gamma(x))) (I guess that's where the definition comes from). I
> find this definition a bit strange, that's not the one I have seen where
> I see it used, but I certainly don't claim to use what would be
> considered as a reference for this stuff (I mostly use log gamma to deal
> with precision problem of gamma in some statistics computation).
>
> cheers,
>
> David
> _______________________________________________
> Scipy-dev mailing list
> Scipy-dev@scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-dev
>
More information about the Scipy-dev
mailing list