Beta function...

[Ben Abbott]

On Monday, March 17, 2008, at 03:45PM, "Ben Abbott" <bpabbott at mac.com> wrote:
>On Monday, March 17, 2008, at 02:56PM, "Jorge Londoño" <jmlon at yahoo.com> wrote:
>>I just happen to jump into a problem that require me to compute the Beta function with negative arguments.
>>If you look at the definition of the beta function, for example here:
>>
>>http://mathworld.wolfram.com/BetaFunction.html, or
>>http://en.wikipedia.org/wiki/Beta_function
>>
>>you will notice it is ok to evaluate the beta function at negative values, as far as they are not integers.
>>
>>The current implementation of the Beta function uses the logarithmic gamma function:
>>
>>retval = exp (gammaln (a) + gammaln (b) - gammaln (a+b))
>>
>>which does not work for negative numbers. (This line is taken from /usr/share/octave/2.9.12/m/specfun/beta.m)
>>
>>Using the alternative definition of the gamma function:
>>x = gamma(a)*gamma(b)/gamma(a+b)
>>
>>it returns the correct values for negative numbers.
>>
>>Jorge Londono
>
>I did a quick comparison of gamma/gammaln
>
>octave:30> x = 1:1e6;
>octave:31> tic;gammaln(x);toc
>Elapsed time is 0.18 seconds.
>octave:32> tic;gamma(x);toc
>Elapsed time is 0.05 seconds.
>
>If speed is not the issue, does anyone have an idea why gammaln was used instead of gamma?
>
>Ben
>

I've answered my own question. The gammaln approach has superior numerical stability.

I also noticed that beta(a,b) returns value even when a and/or b are negative. Although they are not necessarily correct. This is because Octave's gammaln (a) = log (abs (gamma (a))).

Matlab uses an algorithm for gammaln that returns Infs for negative arguments.

For example,

octave:95> gammaln(-.5)
ans =  1.2655
octave:96> log(gamma(-.5))
ans =  1.2655 + 3.1416i

The gammaln result is wrong due to its implementation that apparently implements as log(gamma(abs(x))).

For the same example, Matlab returns

>> gammaln(-.5)
ans =  Inf
>> log(gamma(-.5))
ans = 1.2655 +     3.1416i

Matlab's result, for gammaln, doesn't make any sense to me at all.

Another example using beta

octave:99> beta(-.25,-.25)
ans =  6.7777
octave:100> gamma(-.25)*gamma(-.25)/gamma(-2*0.25)
ans = -6.7777

Matab returns

>> beta(-.25,-.25)
ans =  NaN
>> gamma(-.25)*gamma(-.25)/gamma(-2*0.25)
ans =   -6.7777

I can make a trivial change to beta.m, and replace gammaln(x) with log(gamma(x)). That will resolve the problems for negative arguments. 

But what of gammaln?  Should it be modified as well to return the correct answers for negative inputs?

Anyone?

Ben