tdist.c

This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY without ev

/* cdf/tdist.c
 *
 * Copyright (C) 2002 Jason H. Stover.
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307, USA.
 */

/*
 * Computes the Student's t cumulative distribution function using
 * the method detailed in 
 * 
 * W.J. Kennedy and J.E. Gentle, "Statistical Computing." 1980. 
 * Marcel Dekker. ISBN 0-8247-6898-1.
 *
 * G.W. Hill and A.W. Davis. "Generalized asymptotic expansions
 * of Cornish-Fisher type." Annals of Mathematical Statistics, 
 * vol. 39, 1264-1273. 1968.
 *
 * G.W. Hill. "Algorithm 395: Student's t-distribution," Communications
 * of the ACM, volume 13, number 10, page 617. October 1970.
 *
 * G.W. Hill, "Remark on algorithm 395: Student's t-distribution,"
 * Transactions on Mathematical Software, volume 7, number 2, page
 * 247. June 1982.
 */

#include <config.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>

#include "beta_inc.c"

static double
poly_eval (const double c[], unsigned int n, double x)
{
  unsigned int i;
  double y = c[0] * x;

  for (i = 1; i < n; i++)
    {
      y = x * (y + c[i]);
    }

  y += c[n];

  return y;
}

/* 
 * Use the Cornish-Fisher asymptotic expansion to find a point u such
 * that gsl_cdf_gauss(y) = tcdf(t).
 * 
 */

static double
cornish_fisher (double t, double n)
{
  const double coeffs6[10] = {
    0.265974025974025974026,
    5.449696969696969696970,
    122.20294372294372294372,
    2354.7298701298701298701,
    37625.00902597402597403,
    486996.1392857142857143,
    4960870.65,
    37978595.55,
    201505390.875,
    622437908.625
  };
  const double coeffs5[8] = {
    0.2742857142857142857142,
    4.499047619047619047619,
    78.45142857142857142857,
    1118.710714285714285714,
    12387.6,
    101024.55,
    559494.0,
    1764959.625
  };
  const double coeffs4[6] = {
    0.3047619047619047619048,
    3.752380952380952380952,
    46.67142857142857142857,
    427.5,
    2587.5,
    8518.5
  };
  const double coeffs3[4] = {
    0.4,
    3.3,
    24.0,
    85.5
  };

  double a = n - 0.5;
  double b = 48.0 * a * a;

  double z2 = a * log1p (t * t / n);
  double z = sqrt (z2);

  double p5 = z * poly_eval (coeffs6, 9, z2);
  double p4 = -z * poly_eval (coeffs5, 7, z2);
  double p3 = z * poly_eval (coeffs4, 5, z2);
  double p2 = -z * poly_eval (coeffs3, 3, z2);
  double p1 = z * (z2 + 3.0);
  double p0 = z;

  double y = p5;
  y = (y / b) + p4;
  y = (y / b) + p3;
  y = (y / b) + p2;
  y = (y / b) + p1;
  y = (y / b) + p0;

  if (t < 0)
    y *= -1;

  return y;
}

#if 0
/*
 * Series approximation for t > 4.0. This needs to be fixed;
 * it shouldn't subtract the result from 1.0. A better way is
 * to use two different series expansions. Figuring this out
 * means rummaging through Fisher's paper in Metron, v5, 1926, 
 * "Expansion of Student's integral in powers of n^{-1}."
 */

#define MAXI 40

static double
normal_approx (const double x, const double nu)
{
  double y;
  double num;
  double diff;
  double q;
  int i;
  double lg1, lg2;

  y = 1 / sqrt (1 + x * x / nu);
  num = 1.0;
  q = 0.0;
  diff = 2 * GSL_DBL_EPSILON;
  for (i = 2; (i < MAXI) && (diff > GSL_DBL_EPSILON); i += 2)
    {
      diff = q;
      num *= y * y * (i - 1) / i;
      q += num / (nu + i);
      diff = q - diff;
    }
  q += 1 / nu;
  lg1 = gsl_sf_lngamma (nu / 2.0);
  lg2 = gsl_sf_lngamma ((nu + 1.0) / 2.0);

  diff = lg2 - lg1;
  q *= pow (y, nu) * exp (diff) / sqrt (M_PI);

  return q;
}
#endif

double
gsl_cdf_tdist_P (const double x, const double nu)
{
  double P;

  double x2 = x * x;

  if (nu > 30 && x2 < 10 * nu)
    {
      double u = cornish_fisher (x, nu);
      P = gsl_cdf_ugaussian_P (u);

      return P;
    }

  if (x2 < nu)
    {
      double u = x2 / nu;
      double eps = u / (1 + u);

      if (x >= 0)
        {
          P = beta_inc_AXPY (0.5, 0.5, 0.5, nu / 2.0, eps);
        }
      else
        {
          P = beta_inc_AXPY (-0.5, 0.5, 0.5, nu / 2.0, eps);
        }
    }
  else
    {
      double v = nu / (x * x);
      double eps = v / (1 + v);

      if (x >= 0)
        {
          P = beta_inc_AXPY (-0.5, 1.0, nu / 2.0, 0.5, eps);
        }
      else
        {
          P = beta_inc_AXPY (0.5, 0.0, nu / 2.0, 0.5, eps);
        }
    }

  return P;
}


double
gsl_cdf_tdist_Q (const double x, const double nu)
{
  double Q;

  double x2 = x * x;

  if (nu > 30 && x2 < 10 * nu)
    {
      double u = cornish_fisher (x, nu);
      Q = gsl_cdf_ugaussian_Q (u);

      return Q;
    }

  if (x2 < nu)
    {
      double u = x2 / nu;
      double eps = u / (1 + u);

      if (x >= 0)
        {
          Q = beta_inc_AXPY (-0.5, 0.5, 0.5, nu / 2.0, eps);
        }
      else
        {
          Q = beta_inc_AXPY (0.5, 0.5, 0.5, nu / 2.0, eps);
        }
    }
  else
    {
      double v = nu / (x * x);
      double eps = v / (1 + v);

      if (x >= 0)
        {
          Q = beta_inc_AXPY (0.5, 0.0, nu / 2.0, 0.5, eps);
        }
      else
        {
          Q = beta_inc_AXPY (-0.5, 1.0, nu / 2.0, 0.5, eps);
        }
    }

  return Q;
}