normal.c

it is regression Algorithm in C/C++.

        *x +     5.4637849111641143699)
        *x +     6.6579046435011037772)
      / (((((((    2.04426310338993978564e-15
        *x +     1.4215117583164458887e-7)
        *x +     1.8463183175100546818e-5)
        *x +     7.868691311456132591e-4)
        *x +     0.0148753612908506148525)
        *x +     0.13692988092273580531)
        *x +     0.59983220655588793769)
        *x +     1.0);          /* evaluate the rational function */
  }
  return (prob < 0.5) ? -x : x; /* retransform to right tail if nec. */
}  /* unitqtlP() */

#endif
/*----------------------------------------------------------------------
References: - R.E. Odeh and J.O. Evans.
              The percentage points of the normal distribution.
              Applied Statistics 22:96--97, 1974
            - J.D. Beasley and S.G. Springer.
              Algorithm AS 111:
              The percentage points of the normal distribution.
              Applied Statistics 26:118--121, 1977
            - M.J. Wichura
              Algorithm AS 241:
              The percentage points of the normal distribution.
              Applied Statistics 37:477--484, 1988
----------------------------------------------------------------------*/

double unitrand (double randfn(void))
{                               /* --- compute N(0,1) distrib. number */
  static double b;              /* buffer for random number */
  double x, y, r;               /* coordinates and radius */

  if (b != 0.0) {               /* if the buffer is full, */
    x = b; b = 0; return x; }   /* return the buffered number */
  do {                          /* pick a random point */
    x = 2.0*randfn()-1.0;       /* in the unit square [-1,1]^2 */
    y = 2.0*randfn()-1.0;       /* and check whether it lies */
    r = x*x +y*y;               /* inside the unit circle */
  } while ((r > 1) || (r == 0));
  r = sqrt(-2*log(r)/r);        /* factor for Box-Muller transform */
  b = x *r;                     /* save one of the random numbers */
  return y *r;                  /* and return the other */
}  /* unitrand() */

/*----------------------------------------------------------------------
Source for the polar method to generate normally distributed numbers:
D.E. Knuth.
The Art of Computer Programming, Vol. 2: Seminumerical Algorithms
Addison-Wesley, Reading, MA, USA 1998
pp. 122-123
----------------------------------------------------------------------*/

double normpdf (double x, double mean, double var)
{                               /* --- probability density function */
  assert(var >= 0);             /* check the function arguments */
  if (var == 1) return unitpdf(x-mean);
  return unitpdf((x-mean) /sqrt(var));
}  /* normpdf() */

/*--------------------------------------------------------------------*/

double normcdfP (double x, double mean, double var)
{                               /* --- cumulative distribution fn. */
  assert(var >= 0);             /* check the function arguments */
  if (var == 1) return unitcdf(x-mean);
  return unitcdf((x-mean) /sqrt(var));
}  /* normcdfP() */

/*--------------------------------------------------------------------*/

double normcdfQ (double x, double mean, double var)
{                               /* --- cumulative distribution fn. */
  assert(var >= 0);             /* check the function arguments */
  if (var == 1) return unitcdfP(mean-x);
  return unitcdfP((mean-x) /sqrt(var));
}  /* normcdfQ() */

/*--------------------------------------------------------------------*/

double normqtlP (double prob, double mean, double var)
{                               /* --- quantile of normal distrib. */
  assert((var  >  0)            /* check the function arguments */
      && (prob >= 0) && (prob <= 1));
  if (var == 1) return mean +unitqtl(prob);
  return mean +unitqtlP(prob) *sqrt(var);
}  /* normqtlP() */

/*--------------------------------------------------------------------*/

double normqtlQ (double prob, double mean, double var)
{                               /* --- quantile of normal distrib. */
  assert((var  >  0)            /* check the function arguments */
      && (prob >= 0) && (prob <= 1));
  if (var == 1) return mean -unitqtlP(prob);
  return mean -unitqtlP(prob) *sqrt(var);
}  /* normqtlQ() */

/*--------------------------------------------------------------------*/

double normrand (double randfn(void), double mean, double var)
{                               /* --- cumulative distribution fn. */
  assert(var >= 0);             /* check the function arguments */
  if (var == 1) return mean +unitrand(randfn);
  return mean +unitrand(randfn) *sqrt(var);
}  /* normrand() */

/*----------------------------------------------------------------------
  Main Functions
----------------------------------------------------------------------*/
#ifdef NORMPDF_MAIN

int main (int argc, char *argv[])
{                               /* --- main function */
  double mean = 0;              /* mean value */
  double var  = 1;              /* variance */
  double x;                     /* argument value */

  if ((argc < 2) || (argc > 4)){/* if wrong number of arguments */
    printf("usage: %s arg [mean variance]\n", argv[0]);
    printf("compute probability density function "
           "of the normal distribution\n");
    return 0;                   /* print a usage message */
  }                             /* and abort the program */
  x = atof(argv[1]);            /* get the argument value */
  if (argc > 2) mean = atof(argv[2]);
  if (argc > 3) var  = atof(argv[3]);
  if (var  < 0) {               /* get the parameters */
    printf("%s: invalid variance\n", argv[0]); return -1; }
  printf("normal: f(%.16g;%.16g,%.16g) = %.16g\n",
         x, mean, var, normpdf(x, mean, var));
  return 0;                     /* compute and print density */
}  /* main() */

#endif
/*--------------------------------------------------------------------*/
#ifdef NORMCDF_MAIN

int main (int argc, char *argv[])
{                               /* --- main function */
  double mean = 0;              /* mean value */
  double var  = 1;              /* variance */
  double x;                     /* argument value */

  if ((argc < 2) || (argc > 4)){/* if wrong number of arguments */
    printf("usage: %s arg [mean variance]\n", argv[0]);
    printf("compute cumulative distribution function "
           "of the normal distribution\n");
    return 0;                   /* print a usage message */
  }                             /* and abort the program */
  x = atof(argv[1]);            /* get the argument value */
  if (argc > 2) mean = atof(argv[2]);
  if (argc > 3) var  = atof(argv[3]);
  if (var  < 0) {               /* get the parameters */
    printf("%s: invalid variance\n", argv[0]); return -1; }
  printf("normal: F(% .16g;%.16g,%.16g) = %.16g\n",
         x, mean, var, normcdfP(x, mean, var));
  printf("    1 - F(% .16g;%.16g,%.16g) = %.16g\n",
         x, mean, var, normcdfQ(x, mean, var));
  return 0;                     /* compute and print probability */
}  /* main() */

#endif
/*--------------------------------------------------------------------*/
#ifdef NORMQTL_MAIN

int main (int argc, char *argv[])
{                               /* --- main function */
  double mean = 0;              /* mean value */
  double var  = 1;              /* variance */
  double prob;                  /* probability */

  if ((argc < 2) || (argc > 4)){/* if wrong number of arguments */
    printf("usage: %s prob [mean variance]\n", argv[0]);
    printf("compute quantile of the normal distribution\n");
    return 0;                   /* print a usage message */
  }                             /* and abort the program */
  prob = atof(argv[1]);         /* get the probability */
  if ((prob < 0) || (prob > 1)) {
    printf("%s: invalid probability\n", argv[0]); return -1; }
  if (argc > 2) mean = atof(argv[2]);
  if (argc > 3) var  = atof(argv[3]);
  if (var  < 0) {               /* get the parameters */
    printf("%s: invalid variance\n", argv[0]); return -1; }
  printf("normal: F(% .16g;%.16g,%.16g) = %.16g\n",
         normqtlP(prob, mean, var), mean, var, prob);
  printf("    1 - F(% .16g;%.16g,%.16g) = %.16g\n",
         normqtlQ(prob, mean, var), mean, var, prob);
  return 0;                     /* compute and print quantile */
}  /* main() */

#endif