GRASS GIS 8 Programmer's Manual 8.2.0(2022)-exported
normp.c
Go to the documentation of this file.
1#include <math.h>
2
3
4/*-
5 * SUBROUTINE Cdhc_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
23double Cdhc_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}
double Cdhc_normp(double z)
Definition: normp.c:23