normal.c

it is regression Algorithm in C/C++.

/*----------------------------------------------------------------------
  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 0

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 */
  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)