GRASS GIS 8 Programmer's Manual 8.4.1(2025)-45ca3179ab
Loading...
Searching...
No Matches
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
22double 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}
double Cdhc_normp(double z)
Definition normp.c:22