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