GRASS Programmer's Manual  6.5.svn(2014)-r66266
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
normp.c
Go to the documentation of this file.
1 #include <math.h>
2 
3 
4 /*-
5  * SUBROUTINE NORMP(Z, P, Q, PDF)
6  *
7  * Normal distribution probabilities accurate to 1.e-15.
8  * Z = no. of standard deviations from the mean.
9  * P, Q = probabilities to the left & right of Z. P + Q = 1.
10  * PDF = the probability density.
11  *
12  * Based upon algorithm 5666 for the error function, from:
13  * Hart, J.F. et al, 'Computer Approximations', Wiley 1968
14  *
15  * Programmer: Alan Miller
16  *
17  * Latest revision - 30 March 1986
18  *
19  */
20 
21 /* Conversion to C by James Darrell McCauley, 24 September 1994 */
22 
23 double normp(double z)
24 {
25  static double p[7] = { 220.2068679123761, 221.2135961699311,
26  112.079291497870, 33.91286607838300, 6.37396220353165,
27  0.7003830644436881, 0.352624965998910e-1
28  };
29  static double q[8] = { 440.4137358247522, 793.8265125199484,
30  637.3336333788311, 296.5642487796737, 86.78073220294608,
31  16.06417757920695, 1.755667163182642, 0.8838834764831844e-1
32  };
33  static double cutoff = 7.071, root2pi = 2.506628274631001;
34  double zabs, expntl;
35  double pee, queue, pdf;
36 
37  zabs = fabs(z);
38 
39  if (zabs > 37.0) {
40  pdf = 0.0;
41  if (z > 0.0) {
42  pee = 1.0;
43  queue = 0.0;
44  }
45  else {
46  pee = 0.0;
47  queue = 1.0;
48  }
49 
50  return pee;
51  }
52 
53  expntl = exp(-0.5 * zabs * zabs);
54  pdf = expntl / root2pi;
55 
56  if (zabs < cutoff)
57  pee = expntl * ((((((p[6] * zabs + p[5]) * zabs + p[4])
58  * zabs + p[3]) * zabs + p[2]) * zabs +
59  p[1]) * zabs + p[0])
60  / (((((((q[7] * zabs + q[6]) * zabs + q[5]) * zabs + q[4])
61  * zabs + q[3]) * zabs + q[2]) * zabs + q[1]) * zabs + q[0]);
62  else
63  pee = pdf / (zabs + 1.0 / (zabs + 2.0 / (zabs + 3.0 / (zabs + 4.0 /
64  (zabs +
65  0.65)))));
66 
67  if (z < 0.0)
68  queue = 1.0 - pee;
69  else {
70  queue = pee;
71  pee = 1.0 - queue;
72  }
73 
74  return pee;
75 }
tuple q
Definition: forms.py:2019
double normp(double)
Definition: normp.c:23