GRASS Programmer's Manual  6.5.svn(2014)-r66266
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
as66.c
Go to the documentation of this file.
1 
2 /*-Algorithm AS 66
3  * The Normal Integral, by I. D. Hill, 1973.
4  * Applied Statistics 22(3):424-427.
5  *
6  * Translation to C by James Darrell McCauley, mccauley@ecn.purdue.edu.
7  *
8  * Calculates the upper or lower tail area of the standardized normal
9  * curve corresponding to any given argument.
10  *
11  * x - the argument value
12  * upper: 1 -> the area from x to \infty
13  * 0 -> the area from -\infty to x
14  *
15  * Notes:
16  * The constant LTONE should be set to the value at which the
17  * lower tail area becomes 1.0 to the accuracy of the machine.
18  * LTONE=(n+9)/3 gives the required value accurately enough, for a
19  * machine that produces n decimal digits in its real numbers.
20  *
21  * The constant UTZERO should be set to the value at which the upper
22  * tail area becomes 0.0 to the accuracy of the machine. This may be
23  * taken as the value such that exp(-0.5 * UTZERO * UTZERO) /
24  * (UTZERO * sqrt(2*M_PI)) is just greater than the smallest allowable
25  * real numbers.
26  */
27 
28 #include <math.h>
29 
30 
31 #define LTONE 7.0
32 #define UTZERO 18.66
33 
34 
35 double alnorm(double x, int upper)
36 {
37  double ret, z, y;
38  int up;
39 
40  up = upper;
41  z = x;
42 
43  if (x < 0.0) {
44  up = up == 0 ? 1 : 0;
45  z = -x;
46  }
47 
48  if (!(z <= LTONE || (up == 1 && z <= UTZERO)))
49  ret = 0.0;
50  else {
51  y = 0.5 * z * z;
52  if (z <= 1.28)
53  ret = 0.5 - z * (0.398942280444 - 0.399903438504 * y /
54  (y + 5.75885480458 - 29.8213557808 /
55  (y + 2.62433121679 + 48.6959930692 /
56  (y + 5.92885724438))));
57  else
58  ret = 0.398942280385 * exp(-y) /
59  (z - 3.8052e-8 + 1.00000615302 /
60  (z + 3.98064794e-4 + 1.98615381364 /
61  (z - 0.151679116635 + 5.29330324926 /
62  (z + 4.8385912808 - 15.1508972451 /
63  (z + 0.742380924027 + 30.789933034 /
64  (z + 3.99019417011))))));
65  }
66 
67  if (up == 0)
68  ret = 1.0 - ret;
69 
70  return ret;
71 }
double alnorm(double x, int upper)
Definition: as66.c:35
int y
Definition: plot.c:34
#define LTONE
Definition: as66.c:31
#define UTZERO
Definition: as66.c:32