/*----------------------------------------------------------------------  File    : normal.c  Contents: normal distribution functions  Author  : Christian Borgelt  History : 2002.11.04 file created            2008.03.14 quantile function, main function added            2008.03.15 cumulative distribution function added----------------------------------------------------------------------*/#if defined(NORMPDF_MAIN) \ || defined(NORMCDF_MAIN) \ || defined(NORMQTL_MAIN)#include <stdio.h>#include <stdlib.h>#endif#include <float.h>#include <assert.h>#include <math.h>#include "normal.h"/*----------------------------------------------------------------------  Preprocessor Definitions----------------------------------------------------------------------*/#define SQRT_2      1.41421356237309504880168872421  /* \sqrt(2) */#define _1_SQRT_2PI 0.39894228040143267793994605993  /* 1/\sqrt(2\pi) *//*----------------------------------------------------------------------  Functions----------------------------------------------------------------------*/double unitpdf (double x){                               /* --- probability density function */  return exp(-0.5*x*x) *_1_SQRT_2PI;}  /* unitpdf() *//*--------------------------------------------------------------------*/double unitcdfP (double x){                               /* --- cumulative distribution fn. */  double y, z, u;               /* square, absolute value of x */  if (x >   8.572) return 1.0;  /* if outside proper interval, */  if (x < -37.519) return 0.0;  /* return the limiting values */  z = fabs(x);                  /* exploit the symmetry */  if (z < 0.5*2.2204460492503131e-16)    return 0.5;                 /* treat 0 as a special case */  if (z < 0.66291) {            /* if fairly close to zero */    y = x*x;                   /* compute the square of x */    z = ((((     0.065682337918207449113        *y +     2.2352520354606839287)        *y +   161.02823106855587881)        *y +  1067.6894854603709582)        *y + 18154.981253343561249)      / ((((         y +    47.20258190468824187)        *y +   976.09855173777669322)        *y + 10260.932208618978205)        *y + 45507.789335026729956);    return 0.5 +x *z;           /* compute with rational function */  }                               if (z < 4*SQRT_2) {           /* if medium value */    u = (((((((( 1.0765576773720192317e-8        *z +     0.39894151208813466764)        *z +     8.8831497943883759412)        *z +    93.506656132177855979)        *z +   597.27027639480026226)        *z +  2494.5375852903726711)        *z +  6848.1904505362823326)        *z + 11602.651437647350124)        *z +  9842.7148383839780218)      / ((((((((         z +    22.266688044328115691)        *z +   235.38790178262499861)        *z +  1519.377599407554805)        *z +  6485.558298266760755)        *z + 18615.571640885098091)        *z + 34900.952721145977266)        *z + 38912.003286093271411)        *z + 19685.429676859990727);    z = u *exp(-0.5*x*x); }     /* compute with rational function */  else {                        /* if tail value */    y = 1 /(x*x);               /* compute the inverse square of x */    u = (((((    0.02307344176494017303        *y +     0.21589853405795699)        *y +     0.1274011611602473639)        *y +     0.022235277870649807)        *y +     0.001421619193227893466)        *y +     2.9112874951168792e-5)      / (((((         y +     1.28426009614491121)        *y +     0.468238212480865118)        *y +     0.0659881378689285515)        *y +     0.00378239633202758244)        *y +     7.29751555083966205e-5);    z = (((_1_SQRT_2PI) -y*u) /z) *exp(-0.5*x*x);  }                             /* compute with rational function */  return (x > 0) ? 1-z : z;     /* exploit the symmetry */}  /* unitcdfP() *//*----------------------------------------------------------------------References: - W.J. Cody.              Rational Chebyshev Approximations for the Error Function.              Mathematics of Computation, 23(107):631-637, 1969            - W. Fraser and J.F Hart.              On the Computation of Rational Approximations              to Continuous Functions.              Communications of the ACM 5, 1962            - W.J. Kennedy Jr. and J.E. Gentle.              Statistical Computing.              Marcel Dekker, 1980----------------------------------------------------------------------*/#if 0double unitqtlP (double prob){                               /* --- quantile of normal distrib. */  double p, x;                  /*     with mean 0 and variance 1 */  assert((prob >= 0) && (prob <= 1));  /* check the function argument */  if (prob >= 1.0) return  DBL_MAX; /* check for limiting values */  if (prob <= 0.0) return -DBL_MAX; /* and return extrema */  if (prob == 0.5) return 0;    /* treat prob = 0.5 as a special case */  p = (prob > 0.5) ? 1-prob : prob; /* transform to left tail if nec. */  if (p < 1e-20) return (prob < 0.5) ? -10 : 10;  x = sqrt(log(1/(p*p)));       /* compute quotient of polynomials */  x += ((((-0.453642210148e-4   /* (rational function approximation) */       *x  -0.0204231210245)       *x  -0.342242088547)       *x  -1.0)       *x  -0.322232431088)     / (((( 0.0038560700634       *x  +0.103537752850)       *x  +0.531103462366)       *x  +0.588581570495)       *x  +0.0993484626060);  return (prob < 0.5) ? -x : x; /* retransform to right tail if nec. */}  /* unitqtlP() */#else /*--------------------------------------------------------------*/double unitqtlP (double prob){                               /* --- quantile of normal distrib. */  double p, x;                  /*     with mean 0 and variance 1 */  assert((prob >= 0) && (prob <= 1));  /* check the function argument */  if (prob >= 1.0) return  DBL_MAX; /* check for limiting values */  if (prob <= 0.0) return -DBL_MAX; /* and return extrema */  p = prob -0.5;  if (fabs(p) <= 0.425) {       /* if not tail */    x = 0.180625 - p*p;         /* get argument of rational function */    x = (((((((2509.0809287301226727        *x +  33430.575583588128105)        *x +  67265.770927008700853)        *x +  45921.953931549871457)        *x +  13731.693765509461125)        *x +   1971.5909503065514427)        *x +    133.14166789178437745)        *x +      3.387132872796366608)      / (((((((5226.495278852854561        *x +  28729.085735721942674)        *x +  39307.89580009271061)        *x +  21213.794301586595867)        *x +   5394.1960214247511077)        *x +    687.1870074920579083)        *x +     42.313330701600911252)        *x +      1.0);         /* evaluate the rational function */    return p *x;                /* and return the computed value */  }  p = (prob > 0.5) ? 1-prob : prob;  x = sqrt(-log(p));            /* transform to left tail if nec. */  if (x <= 5) {                 /* if not extreme tail */    x -= 1.6;                   /* get argument of rational function */    x = (((((((  7.7454501427834140764e-4        *x +     0.0227238449892691845833)        *x +     0.24178072517745061177)        *x +     1.27045825245236838258)        *x +     3.64784832476320460504)        *x +     5.7694972214606914055)        *x +     4.6303378461565452959)        *x +     1.42343711074968357734)      / (((((((  1.05075007164441684324e-9        *x +     5.475938084995344946e-4)        *x +     0.0151986665636164571966)        *x +     0.14810397642748007459)        *x +     0.68976733498510000455)        *x +     1.6763848301838038494)        *x +     2.05319162663775882187)        *x +     1.0); }        /* evaluate the rational function */  else {                        /* if extreme tail */    x -= 5;                     /* get argument of rational function */    x = (((((((  2.01033439929228813265e-7        *x +     2.71155556874348757815e-5)        *x +     0.0012426609473880784386)        *x +     0.026532189526576123093)        *x +     0.29656057182850489123)        *x +     1.7848265399172913358)