Re: Generalised Gamma Distribution

[Martin Jansche <jansche at ling dot ohio-state dot edu>]

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