Generalised Gamma Distribution

Generalised Gamma Distribution

[Adam Johansen]

Hello,

Am I correct in thinking that the GSL doesn't contain an implementation of
the generalised gamma function and that there is no implementation of
the gamma function for complex variables?

If so, does anybody have any advice on implementing one?

Thanks in advance,
Adam Johansen

Re: Generalised Gamma Distribution

[Brian Gough]

Adam Johansen writes:
 > Hello,
 > 
 > Am I correct in thinking that the GSL doesn't contain an
 > implementation of the generalised gamma function

I'm not familiar with the definition, which function does it
correspond to in Abramowitz and Stegun?

 > and that there is no implementation of the gamma function for
 > complex variables?

There is gsl_sf_lngamma_complex_e for Log[Gamma(z)].

 > If so, does anybody have any advice on implementing one?

There is some general advice in the design document at
http://sources.redhat.com/gsl/

-- 
Brian Gough

Re: Generalised Gamma Distribution

[Martin Jansche]

On Tue, 14 Jan 2003, Brian Gough wrote:

> Adam Johansen writes:
>  >
>  > Am I correct in thinking that the GSL doesn't contain an
>  > implementation of the generalised gamma function
>
> I'm not familiar with the definition, which function does it
> correspond to in Abramowitz and Stegun?

I think Adam was referring to 6.5.2 $\gamma(a,x)$ and 6.5.3
$\Gamma(a,x)$ on page 260.  If that's what he meant then I agree with
him that GSL does not contain an implementation of either 6.5.2 or
6.5.3 (the incomplete (G/g)amma functions), only of 6.5.1 $P(a,x)$
(the regularized/normalized incomplete gamma function) and $Q(a,x) =
1 - P(a,x)$.

I don't think that one can generally obtain Gamma(a,x) from Q(a,x) and
Gamma(a) because of singularities.  For example, Gamma(0,1) =
expint_E1(1) ~= 0.21938393; or Gamma(-1, 1) ~= 0.14849551.

If you need Gamma(a,x), here's a poor substitute with all sorts of
silly aspects (like returning Re(Gamma(a,x)) for x<0, a<0 and
a an integer, but returning NaN if x<0, a<0 and a is not an integer):

double
poor_persons_Gamma(double a, double x)
{
  double aint = floor(a + 0.5);
  if ((fabs(a - aint) < MY_EPS)
      && (aint >= INT_MIN) && (aint <= INT_MAX))
    /* a is an integer */
    {
      int p = (int) aint;
      if (p == 1)
	return exp (-x);
      else if (p == 0)
	return gsl_sf_expint_E1 (x);
      else if (p < 0)
	{
	  /* A&S 6.5.19 */
	  int n = -p;
	  double g = 0.0;
	  int j;
	  for (j = 0; j < n; j++)
	    g += ((j&1)? -1 : 1) * gsl_sf_fact (j) / pow (x, j+1);
	  g *= exp (-x);
	  g = gsl_sf_expint_E1 (x) - g;
	  g /= gsl_sf_fact (n);
	  g *= (j&1)? -1 : 1;
	  return g;
	}
    }
  if (fabs(x) < MY_EPS)
    /* x ~= 0 */
    {
      return gsl_sf_gamma (a);
    }
  else
    {
      double g = gsl_sf_hyperg_1F1 (a, a+1, -x);
      g *= - pow (x, a) / a;
      g += gsl_sf_gamma(a);
      return g;
    }
}

Somebody ought to do this right.

BTW the incomplete Beta function is missing too, and the situation is
quite similar.

- martin

Re: Generalised Gamma Function

[Adam Johansen]

Many Thanks for your suggestions.

I apologies for not being clearer about what I was referring to -- and for
typing "distribution" whilst thinking "function" in the title.

> On Tue, 14 Jan 2003, Brian Gough wrote:
> > I'm not familiar with the definition, which function does it
> > correspond to in Abramowitz and Stegun?

I'm afraid that I don't thyink that there is a definition for it in
Abramowitz and Stegun. In actual fact it can be obtained trivially from
the Gamma function. The definition in Bernardo and Smith, "Bayesian
Theory", Wiley, 2002 p138 is:

Gamma_k(a) = Pi^(k(k-1)/4) * product(i = 1..k, Gamma(0.5 * (2a + 1 - i)))

Sorry if this is an unusual name for this function, or you know it as
something different.

I hope that's intelligible. This is, of course, a trivial product of
normal Gamma functions if the parameter a is restricted to real numbers.
I've implemented this, but I assume it's of no interest to anyone. My
application does, it transpires, only really require this case so I'm not
going to complicate things any further at this stage.

Brian: thanks for pointing out the gsl_sf_lngamma_complex_e function. I'm
not sure how I missed it, but this was, in fact pretty much exactly what I
was looking for.

Apologies for any confusion.

Adam