gamma.c

it is regression Algorithm in C/C++.

  assert((n > 0) && (x > 0));   /* check the function arguments */
  return _series(n, x) *exp(n *log(x) -x);
}  /* lowerGamma() */

/*--------------------------------------------------------------------*/

double upperGamma (double n, double x)
{                               /* --- upper incomplete Gamma fn. */
  assert((n > 0) && (x > 0));   /* check the function arguments */
  return _cfrac(n, x) *exp(n *log(x) -x);
}  /* upperGamma() */

/*--------------------------------------------------------------------*/

double GammaP (double n, double x)
{                               /* --- regularized Gamma function P */
  assert((n > 0) && (x >= 0));  /* check the function arguments */
  if (x <=  0) return 0;        /* treat x = 0 as a special case */
  if (x < n+1) return _series(n, x) *exp(n *log(x) -x -logGamma(n));
  return 1 -_cfrac(n, x) *exp(n *log(x) -x -logGamma(n));
}  /* GammaP() */

/*--------------------------------------------------------------------*/

double GammaQ (double n, double x)
{                               /* --- regularized Gamma function Q */
  assert((n > 0) && (x >= 0));  /* check the function arguments */
  if (x <=  0) return 1;        /* treat x = 0 as a special case */
  if (x < n+1) return 1 -_series(n, x) *exp(n *log(x) -x -logGamma(n));
  return _cfrac(n, x) *exp(n *log(x) -x -logGamma(n));
}  /* GammaQ() */

/*----------------------------------------------------------------------
P(a,x) is also called the regularized gamma function, Q(a,x) = 1-P(a,x).
P(k/2,x/2), where k is a natural number, is the cumulative distribution
function (cdf) of a chi^2 distribution with k degrees of freedom.
----------------------------------------------------------------------*/

double Gammapdf (double x, double k, double theta)
{                               /* --- probability density function */
  assert((k > 0) && (theta > 0));
  if (x <  0) return 0;         /* support is non-negative x */
  if (x <= 0) return (k == 1) ? 1/theta : 0;
  if (k == 1) return exp(-x/theta) /theta;
  return exp ((k-1) *log(x/theta) -x/theta -logGamma(k)) /theta;
}  /* Gammapdf() */

/*--------------------------------------------------------------------*/
#ifdef GAMMAQTL

double GammaqtlP (double prob, double k, double theta)
{                               /* --- quantile of Gamma distribution */
  int    n = 0;                 /* loop variable */
  double x, f, a, d, dx, dp;    /* buffers */

  assert((k > 0) && (theta > 0) /* check the function arguments */
      && (prob >= 0) && (prob <= 1));
  if (prob >= 1.0) return DBL_MAX;
  if (prob <= 0.0) return 0;    /* handle limiting values */
  if      (prob < 0.05) x = exp(logGamma(k) +log(prob) /k);
  else if (prob > 0.95) x = logGamma(k) -log1p(-prob);
  else {                        /* distinguish three prob. ranges */
    f = unitqtlP(prob); a = sqrt(k);
    x = (f >= -a) ? a *f +k : k;
  }                             /* compute initial approximation */
  do {                          /* Lagrange's interpolation */
    dp = prob -GammacdfP(x, k, 1);
    if ((dp == 0) || (++n > 33)) break;
    f = Gammapdf(x, k, 1);
    a = 2 *fabs(dp/x);
    a = dx = dp /((a > f) ? a : f);
    d = -0.25 *((k-1)/x -1) *a*a;
    if (fabs(d) < fabs(a)) dx += d;
    if (x +dx > 0) x += dx;
    else           x /= 2;
  } while (fabs(a) > 1e-10 *x);
  if (fabs(dp) > EPS_QTL *prob) return -1;
  return x *theta;              /* check for convergence and */
}  /* GammaqtlP() */            /* return the computed quantile */

/*--------------------------------------------------------------------*/

double GammaqtlQ (double prob, double k, double theta)
{                               /* --- quantile of Gamma distribution */
  int    n = 0;                 /* loop variable */
  double x, f, a, d, dx, dp;    /* buffers */

  assert((k > 0) && (theta > 0) /* check the function arguments */
      && (prob >= 0) && (prob <= 1));
  if (prob <= 0.0) return DBL_MAX;
  if (prob >= 1.0) return 0;    /* handle limiting values */
  if      (prob < 0.05) x = logGamma(k) -log(prob);
  else if (prob > 0.95) x = exp(logGamma(k) +log1p(-prob) /k);
  else {                        /* distinguish three prob. ranges */
    f = unitqtlQ(prob); a = sqrt(k);
    x = (f >= -a) ? a *f +k : k;
  }                             /* compute initial approximation */
  do {                          /* Lagrange's interpolation */
    dp = prob -GammacdfQ(x, k, 1);
    if ((dp == 0) || (++n > 33)) break;
    f = Gammapdf(x, k, 1);
    a = 2 *fabs(dp/x);
    a = dx = -dp /((a > f) ? a : f);
    d = -0.25 *((k-1)/x -1) *a*a;
    if (fabs(d) < fabs(a)) dx += d;
    if (x +dx > 0) x += dx;
    else           x /= 2;
  } while (fabs(a) > 1e-10 *x);
  if (fabs(dp) > EPS_QTL *prob) return -1;
  return x *theta;              /* check for convergence and */
}  /* GammaqtlQ() */            /* return the computed quantile */

#endif
/*--------------------------------------------------------------------*/
#ifdef GAMMA_MAIN

int main (int argc, char *argv[])
{                               /* --- main function */
  double x;                     /* argument */

  if (argc != 2) {              /* if wrong number of arguments given */
    printf("usage: %s x\n", argv[0]);
    printf("compute (logarithm of) Gamma function\n");
    return 0;                   /* print a usage message */
  }                             /* and abort the program */
  x = atof(argv[1]);            /* get argument */
  if (x <= 0) { printf("%s: x must be > 0\n", argv[0]); return -1; }
  printf("   Gamma(%.16g)  = % .20g\n", x, Gamma(x));
  printf("ln(Gamma(%.16g)) = % .20g\n", x, logGamma(x));
  return 0;                     /* compute and print Gamma function */
}  /* main() */

#endif
/*--------------------------------------------------------------------*/
#ifdef GAMMAPDF_MAIN

int main (int argc, char *argv[])
{                               /* --- main function */
  double shape = 1;             /* shape parameter */
  double scale = 1;             /* scale parameter */
  double x;                     /* argument value */

  if ((argc < 2) || (argc > 4)){/* if wrong number of arguments */
    printf("usage: %s arg [shape scale]\n", argv[0]);
    printf("compute probability density function "
           "of the gamma distribution\n");
    return 0;                   /* print a usage message */
  }                             /* and abort the program */
  x = atof(argv[1]);            /* get the argument value */
  if (argc > 2) shape = atof(argv[2]);
  if (shape <= 0) {             /* get the parameters */
    printf("%s: invalid shape parameter\n", argv[0]); return -1; }
  if (argc > 3) scale = atof(argv[3]);
  if (scale <= 0) {             /* get the parameters */
    printf("%s: invalid scale parameter\n", argv[0]); return -1; }
  printf("gamma: f(%.16g;%.16g,%.16g) = %.16g\n",
         x, shape, scale, Gammapdf(x, shape, scale));
  return 0;                     /* compute and print density */
}  /* main() */

#endif
/*--------------------------------------------------------------------*/
#ifdef GAMMACDF_MAIN

int main (int argc, char *argv[])
{                               /* --- main function */
  double shape = 1;             /* shape parameter */
  double scale = 1;             /* scale parameter */
  double x;                     /* argument value */

  if ((argc < 2) || (argc > 4)){/* if wrong number of arguments */
    printf("usage: %s arg [shape scale]\n", argv[0]);
    printf("compute cumulative distribution function "
           "of the gamma distribution\n");
    return 0;                   /* print a usage message */
  }                             /* and abort the program */
  x = atof(argv[1]);            /* get the argument value */
  if (argc > 2) shape = atof(argv[2]);
  if (shape <= 0) {             /* get the parameters */
    printf("%s: invalid shape parameter\n", argv[0]); return -1; }
  if (argc > 3) scale = atof(argv[3]);
  if (scale <= 0) {             /* get the parameters */
    printf("%s: invalid scale parameter\n", argv[0]); return -1; }
  printf("gamma: F(% .16g;%.16g,%.16g) = %.16g\n",
         x, shape, scale, GammacdfP(x, shape, scale));
  printf("   1 - F(% .16g;%.16g,%.16g) = %.16g\n",
         x, shape, scale, GammacdfQ(x, shape, scale));
  return 0;                     /* compute and print probability */
}  /* main() */

#endif
/*--------------------------------------------------------------------*/
#ifdef GAMMAQTL_MAIN

int main (int argc, char *argv[])
{                               /* --- main function */
  double shape = 1;             /* shape parameter */
  double scale = 1;             /* scale parameter */
  double prob;                  /* argument value */

  if ((argc < 2) || (argc > 4)){/* if wrong number of arguments */
    printf("usage: %s prob [shape scale]\n", argv[0]);
    printf("compute quantile of the gamma distribution\n");
    return 0;                   /* print a usage message */
  }                             /* and abort the program */
  prob = atof(argv[1]);         /* get the probability */
  if ((prob < 0) || (prob > 1)){/* and check it */
    printf("%s: invalid probability\n", argv[0]); return -1; }
  if (argc > 2) shape = atof(argv[2]);
  if (shape <= 0) {             /* get the parameters */
    printf("%s: invalid shape parameter\n", argv[0]); return -1; }
  if (argc > 3) scale = atof(argv[3]);
  if (scale <= 0) {             /* get the parameters */
    printf("%s: invalid scale parameter\n", argv[0]); return -1; }
  printf("gamma: F(% .16g;%.16g,%.16g) = %.16g\n",
         GammaqtlP(prob, shape, scale), shape, scale, prob);
  printf("   1 - F(% .16g;%.16g,%.16g) = %.16g\n",
         GammaqtlQ(prob, shape, scale), shape, scale, prob);
  return 0;                     /* compute and print probability */
}  /* main() */

#endif