htswanalysis/MACS/lib/gsl/gsl-1.11/cdf/tdist.c

   1 /* cdf/tdist.c
   2  *
   3  * Copyright (C) 2002 Jason H. Stover.
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307, USA.
  18  */
  19
  20 /*
  21  * Computes the Student's t cumulative distribution function using
  22  * the method detailed in
  23  *
  24  * W.J. Kennedy and J.E. Gentle, "Statistical Computing." 1980.
  25  * Marcel Dekker. ISBN 0-8247-6898-1.
  26  *
  27  * G.W. Hill and A.W. Davis. "Generalized asymptotic expansions
  28  * of Cornish-Fisher type." Annals of Mathematical Statistics,
  29  * vol. 39, 1264-1273. 1968.
  30  *
  31  * G.W. Hill. "Algorithm 395: Student's t-distribution," Communications
  32  * of the ACM, volume 13, number 10, page 617. October 1970.
  33  *
  34  * G.W. Hill, "Remark on algorithm 395: Student's t-distribution,"
  35  * Transactions on Mathematical Software, volume 7, number 2, page
  36  * 247. June 1982.
  37  */
  38
  39 #include <config.h>
  40 #include <gsl/gsl_cdf.h>
  41 #include <gsl/gsl_sf_gamma.h>
  42 #include <gsl/gsl_math.h>
  43 #include <gsl/gsl_errno.h>
  44
  45 #include "beta_inc.c"
  46
  47 static double
  48 poly_eval (const double c[], unsigned int n, double x)
  49 {
  50   unsigned int i;
  51   double y = c[0] * x;
  52
  53   for (i = 1; i < n; i++)
  54     {
  55       y = x * (y + c[i]);
  56     }
  57
  58   y += c[n];
  59
  60   return y;
  61 }
  62
  63 /*
  64  * Use the Cornish-Fisher asymptotic expansion to find a point u such
  65  * that gsl_cdf_gauss(y) = tcdf(t).
  66  *
  67  */
  68
  69 static double
  70 cornish_fisher (double t, double n)
  71 {
  72   const double coeffs6[10] = {
  73     0.265974025974025974026,
  74     5.449696969696969696970,
  75     122.20294372294372294372,
  76     2354.7298701298701298701,
  77     37625.00902597402597403,
  78     486996.1392857142857143,
  79     4960870.65,
  80     37978595.55,
  81     201505390.875,
  82     622437908.625
  83   };
  84   const double coeffs5[8] = {
  85     0.2742857142857142857142,
  86     4.499047619047619047619,
  87     78.45142857142857142857,
  88     1118.710714285714285714,
  89     12387.6,
  90     101024.55,
  91     559494.0,
  92     1764959.625
  93   };
  94   const double coeffs4[6] = {
  95     0.3047619047619047619048,
  96     3.752380952380952380952,
  97     46.67142857142857142857,
  98     427.5,
  99     2587.5,
 100     8518.5
 101   };
 102   const double coeffs3[4] = {
 103     0.4,
 104     3.3,
 105     24.0,
 106     85.5
 107   };
 108
 109   double a = n - 0.5;
 110   double b = 48.0 * a * a;
 111
 112   double z2 = a * log1p (t * t / n);
 113   double z = sqrt (z2);
 114
 115   double p5 = z * poly_eval (coeffs6, 9, z2);
 116   double p4 = -z * poly_eval (coeffs5, 7, z2);
 117   double p3 = z * poly_eval (coeffs4, 5, z2);
 118   double p2 = -z * poly_eval (coeffs3, 3, z2);
 119   double p1 = z * (z2 + 3.0);
 120   double p0 = z;
 121
 122   double y = p5;
 123   y = (y / b) + p4;
 124   y = (y / b) + p3;
 125   y = (y / b) + p2;
 126   y = (y / b) + p1;
 127   y = (y / b) + p0;
 128
 129   if (t < 0)
 130     y *= -1;
 131
 132   return y;
 133 }
 134
 135 #if 0
 136 /*
 137  * Series approximation for t > 4.0. This needs to be fixed;
 138  * it shouldn't subtract the result from 1.0. A better way is
 139  * to use two different series expansions. Figuring this out
 140  * means rummaging through Fisher's paper in Metron, v5, 1926,
 141  * "Expansion of Student's integral in powers of n^{-1}."
 142  */
 143
 144 #define MAXI 40
 145
 146 static double
 147 normal_approx (const double x, const double nu)
 148 {
 149   double y;
 150   double num;
 151   double diff;
 152   double q;
 153   int i;
 154   double lg1, lg2;
 155
 156   y = 1 / sqrt (1 + x * x / nu);
 157   num = 1.0;
 158   q = 0.0;
 159   diff = 2 * GSL_DBL_EPSILON;
 160   for (i = 2; (i < MAXI) && (diff > GSL_DBL_EPSILON); i += 2)
 161     {
 162       diff = q;
 163       num *= y * y * (i - 1) / i;
 164       q += num / (nu + i);
 165       diff = q - diff;
 166     }
 167   q += 1 / nu;
 168   lg1 = gsl_sf_lngamma (nu / 2.0);
 169   lg2 = gsl_sf_lngamma ((nu + 1.0) / 2.0);
 170
 171   diff = lg2 - lg1;
 172   q *= pow (y, nu) * exp (diff) / sqrt (M_PI);
 173
 174   return q;
 175 }
 176 #endif
 177
 178 double
 179 gsl_cdf_tdist_P (const double x, const double nu)
 180 {
 181   double P;
 182
 183   double x2 = x * x;
 184
 185   if (nu > 30 && x2 < 10 * nu)
 186     {
 187       double u = cornish_fisher (x, nu);
 188       P = gsl_cdf_ugaussian_P (u);
 189
 190       return P;
 191     }
 192
 193   if (x2 < nu)
 194     {
 195       double u = x2 / nu;
 196       double eps = u / (1 + u);
 197
 198       if (x >= 0)
 199         {
 200           P = beta_inc_AXPY (0.5, 0.5, 0.5, nu / 2.0, eps);
 201         }
 202       else
 203         {
 204           P = beta_inc_AXPY (-0.5, 0.5, 0.5, nu / 2.0, eps);
 205         }
 206     }
 207   else
 208     {
 209       double v = nu / (x * x);
 210       double eps = v / (1 + v);
 211
 212       if (x >= 0)
 213         {
 214           P = beta_inc_AXPY (-0.5, 1.0, nu / 2.0, 0.5, eps);
 215         }
 216       else
 217         {
 218           P = beta_inc_AXPY (0.5, 0.0, nu / 2.0, 0.5, eps);
 219         }
 220     }
 221
 222   return P;
 223 }
 224
 225
 226 double
 227 gsl_cdf_tdist_Q (const double x, const double nu)
 228 {
 229   double Q;
 230
 231   double x2 = x * x;
 232
 233   if (nu > 30 && x2 < 10 * nu)
 234     {
 235       double u = cornish_fisher (x, nu);
 236       Q = gsl_cdf_ugaussian_Q (u);
 237
 238       return Q;
 239     }
 240
 241   if (x2 < nu)
 242     {
 243       double u = x2 / nu;
 244       double eps = u / (1 + u);
 245
 246       if (x >= 0)
 247         {
 248           Q = beta_inc_AXPY (-0.5, 0.5, 0.5, nu / 2.0, eps);
 249         }
 250       else
 251         {
 252           Q = beta_inc_AXPY (0.5, 0.5, 0.5, nu / 2.0, eps);
 253         }
 254     }
 255   else
 256     {
 257       double v = nu / (x * x);
 258       double eps = v / (1 + v);
 259
 260       if (x >= 0)
 261         {
 262           Q = beta_inc_AXPY (0.5, 0.0, nu / 2.0, 0.5, eps);
 263         }
 264       else
 265         {
 266           Q = beta_inc_AXPY (-0.5, 1.0, nu / 2.0, 0.5, eps);
 267         }
 268     }
 269
 270   return Q;
 271 }