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

   1 /* cdf/gauss.c
   2  *
   3  * Copyright (C) 2002, 2004 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 cumulative distribution function for the Gaussian
  22  * distribution using a rational function approximation.  The
  23  * computation is for the standard Normal distribution, i.e., mean 0
  24  * and standard deviation 1. If you want to compute Pr(X < t) for a
  25  * Gaussian random variable X with non-zero mean m and standard
  26  * deviation sd not equal to 1, find gsl_cdf_ugaussian ((t-m)/sd).
  27  * This approximation is accurate to at least double precision. The
  28  * accuracy was verified with a pari-gp script.  The largest error
  29  * found was about 1.4E-20. The coefficients were derived by Cody.
  30  *
  31  * References:
  32  *
  33  * W.J. Cody. "Rational Chebyshev Approximations for the Error
  34  * Function," Mathematics of Computation, v23 n107 1969, 631-637.
  35  *
  36  * W. Fraser, J.F Hart. "On the Computation of Rational Approximations
  37  * to Continuous Functions," Communications of the ACM, v5 1962.
  38  *
  39  * W.J. Kennedy Jr., J.E. Gentle. "Statistical Computing." Marcel Dekker. 1980.
  40  *
  41  *
  42  */
  43
  44 #include <config.h>
  45 #include <math.h>
  46 #include <gsl/gsl_math.h>
  47 #include <gsl/gsl_cdf.h>
  48
  49 #ifndef M_1_SQRT2PI
  50 #define M_1_SQRT2PI (M_2_SQRTPI * M_SQRT1_2 / 2.0)
  51 #endif
  52
  53 #define SQRT32 (4.0 * M_SQRT2)
  54
  55 /*
  56  * IEEE double precision dependent constants.
  57  *
  58  * GAUSS_EPSILON: Smallest positive value such that
  59  *                gsl_cdf_gaussian(x) > 0.5.
  60  * GAUSS_XUPPER: Largest value x such that gsl_cdf_gaussian(x) < 1.0.
  61  * GAUSS_XLOWER: Smallest value x such that gsl_cdf_gaussian(x) > 0.0.
  62  */
  63
  64 #define GAUSS_EPSILON  (GSL_DBL_EPSILON / 2)
  65 #define GAUSS_XUPPER (8.572)
  66 #define GAUSS_XLOWER (-37.519)
  67
  68 #define GAUSS_SCALE (16.0)
  69
  70 static double
  71 get_del (double x, double rational)
  72 {
  73   double xsq = 0.0;
  74   double del = 0.0;
  75   double result = 0.0;
  76
  77   xsq = floor (x * GAUSS_SCALE) / GAUSS_SCALE;
  78   del = (x - xsq) * (x + xsq);
  79   del *= 0.5;
  80
  81   result = exp (-0.5 * xsq * xsq) * exp (-1.0 * del) * rational;
  82
  83   return result;
  84 }
  85
  86 /*
  87  * Normal cdf for fabs(x) < 0.66291
  88  */
  89 static double
  90 gauss_small (const double x)
  91 {
  92   unsigned int i;
  93   double result = 0.0;
  94   double xsq;
  95   double xnum;
  96   double xden;
  97
  98   const double a[5] = {
  99     2.2352520354606839287,
 100     161.02823106855587881,
 101     1067.6894854603709582,
 102     18154.981253343561249,
 103     0.065682337918207449113
 104   };
 105   const double b[4] = {
 106     47.20258190468824187,
 107     976.09855173777669322,
 108     10260.932208618978205,
 109     45507.789335026729956
 110   };
 111
 112   xsq = x * x;
 113   xnum = a[4] * xsq;
 114   xden = xsq;
 115
 116   for (i = 0; i < 3; i++)
 117     {
 118       xnum = (xnum + a[i]) * xsq;
 119       xden = (xden + b[i]) * xsq;
 120     }
 121
 122   result = x * (xnum + a[3]) / (xden + b[3]);
 123
 124   return result;
 125 }
 126
 127 /*
 128  * Normal cdf for 0.66291 < fabs(x) < sqrt(32).
 129  */
 130 static double
 131 gauss_medium (const double x)
 132 {
 133   unsigned int i;
 134   double temp = 0.0;
 135   double result = 0.0;
 136   double xnum;
 137   double xden;
 138   double absx;
 139
 140   const double c[9] = {
 141     0.39894151208813466764,
 142     8.8831497943883759412,
 143     93.506656132177855979,
 144     597.27027639480026226,
 145     2494.5375852903726711,
 146     6848.1904505362823326,
 147     11602.651437647350124,
 148     9842.7148383839780218,
 149     1.0765576773720192317e-8
 150   };
 151   const double d[8] = {
 152     22.266688044328115691,
 153     235.38790178262499861,
 154     1519.377599407554805,
 155     6485.558298266760755,
 156     18615.571640885098091,
 157     34900.952721145977266,
 158     38912.003286093271411,
 159     19685.429676859990727
 160   };
 161
 162   absx = fabs (x);
 163
 164   xnum = c[8] * absx;
 165   xden = absx;
 166
 167   for (i = 0; i < 7; i++)
 168     {
 169       xnum = (xnum + c[i]) * absx;
 170       xden = (xden + d[i]) * absx;
 171     }
 172
 173   temp = (xnum + c[7]) / (xden + d[7]);
 174
 175   result = get_del (x, temp);
 176
 177   return result;
 178 }
 179
 180 /*
 181  * Normal cdf for
 182  * {sqrt(32) < x < GAUSS_XUPPER} union { GAUSS_XLOWER < x < -sqrt(32) }.
 183  */
 184 static double
 185 gauss_large (const double x)
 186 {
 187   int i;
 188   double result;
 189   double xsq;
 190   double temp;
 191   double xnum;
 192   double xden;
 193   double absx;
 194
 195   const double p[6] = {
 196     0.21589853405795699,
 197     0.1274011611602473639,
 198     0.022235277870649807,
 199     0.001421619193227893466,
 200     2.9112874951168792e-5,
 201     0.02307344176494017303
 202   };
 203   const double q[5] = {
 204     1.28426009614491121,
 205     0.468238212480865118,
 206     0.0659881378689285515,
 207     0.00378239633202758244,
 208     7.29751555083966205e-5
 209   };
 210
 211   absx = fabs (x);
 212   xsq = 1.0 / (x * x);
 213   xnum = p[5] * xsq;
 214   xden = xsq;
 215
 216   for (i = 0; i < 4; i++)
 217     {
 218       xnum = (xnum + p[i]) * xsq;
 219       xden = (xden + q[i]) * xsq;
 220     }
 221
 222   temp = xsq * (xnum + p[4]) / (xden + q[4]);
 223   temp = (M_1_SQRT2PI - temp) / absx;
 224
 225   result = get_del (x, temp);
 226
 227   return result;
 228 }
 229
 230 double
 231 gsl_cdf_ugaussian_P (const double x)
 232 {
 233   double result;
 234   double absx = fabs (x);
 235
 236   if (absx < GAUSS_EPSILON)
 237     {
 238       result = 0.5;
 239       return result;
 240     }
 241   else if (absx < 0.66291)
 242     {
 243       result = 0.5 + gauss_small (x);
 244       return result;
 245     }
 246   else if (absx < SQRT32)
 247     {
 248       result = gauss_medium (x);
 249
 250       if (x > 0.0)
 251         {
 252           result = 1.0 - result;
 253         }
 254
 255       return result;
 256     }
 257   else if (x > GAUSS_XUPPER)
 258     {
 259       result = 1.0;
 260       return result;
 261     }
 262   else if (x < GAUSS_XLOWER)
 263     {
 264       result = 0.0;
 265       return result;
 266     }
 267   else
 268     {
 269       result = gauss_large (x);
 270
 271       if (x > 0.0)
 272         {
 273           result = 1.0 - result;
 274         }
 275     }
 276
 277   return result;
 278 }
 279
 280 double
 281 gsl_cdf_ugaussian_Q (const double x)
 282 {
 283   double result;
 284   double absx = fabs (x);
 285
 286   if (absx < GAUSS_EPSILON)
 287     {
 288       result = 0.5;
 289       return result;
 290     }
 291   else if (absx < 0.66291)
 292     {
 293       result = gauss_small (x);
 294
 295       if (x < 0.0)
 296         {
 297           result = fabs (result) + 0.5;
 298         }
 299       else
 300         {
 301           result = 0.5 - result;
 302         }
 303
 304       return result;
 305     }
 306   else if (absx < SQRT32)
 307     {
 308       result = gauss_medium (x);
 309
 310       if (x < 0.0)
 311         {
 312           result = 1.0 - result;
 313         }
 314
 315       return result;
 316     }
 317   else if (x > -(GAUSS_XLOWER))
 318     {
 319       result = 0.0;
 320       return result;
 321     }
 322   else if (x < -(GAUSS_XUPPER))
 323     {
 324       result = 1.0;
 325       return result;
 326     }
 327   else
 328     {
 329       result = gauss_large (x);
 330
 331       if (x < 0.0)
 332         {
 333           result = 1.0 - result;
 334         }
 335
 336     }
 337
 338   return result;
 339 }
 340
 341 double
 342 gsl_cdf_gaussian_P (const double x, const double sigma)
 343 {
 344   return gsl_cdf_ugaussian_P (x / sigma);
 345 }
 346
 347 double
 348 gsl_cdf_gaussian_Q (const double x, const double sigma)
 349 {
 350   return gsl_cdf_ugaussian_Q (x / sigma);
 351 }