This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
Re: Generalised Gamma Distribution
- From: Martin Jansche <jansche at ling dot ohio-state dot edu>
- To: <gsl-discuss at sources dot redhat dot com>
- Date: Wed, 15 Jan 2003 02:26:44 -0500 (EST)
- Subject: Re: Generalised Gamma Distribution
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