GRASS Programmer's Manual
6.4.2(2012)
|
00001 00002 /*-Algorithm AS 66 00003 * The Normal Integral, by I. D. Hill, 1973. 00004 * Applied Statistics 22(3):424-427. 00005 * 00006 * Translation to C by James Darrell McCauley, mccauley@ecn.purdue.edu. 00007 * 00008 * Calculates the upper or lower tail area of the standardized normal 00009 * curve corresponding to any given argument. 00010 * 00011 * x - the argument value 00012 * upper: 1 -> the area from x to \infty 00013 * 0 -> the area from -\infty to x 00014 * 00015 * Notes: 00016 * The constant LTONE should be set to the value at which the 00017 * lower tail area becomes 1.0 to the accuracy of the machine. 00018 * LTONE=(n+9)/3 gives the required value accurately enough, for a 00019 * machine that produces n decimal digits in its real numbers. 00020 * 00021 * The constant UTZERO should be set to the value at which the upper 00022 * tail area becomes 0.0 to the accuracy of the machine. This may be 00023 * taken as the value such that exp(-0.5 * UTZERO * UTZERO) / 00024 * (UTZERO * sqrt(2*M_PI)) is just greater than the smallest allowable 00025 * real numbers. 00026 */ 00027 00028 #include <math.h> 00029 00030 00031 #define LTONE 7.0 00032 #define UTZERO 18.66 00033 00034 00035 double alnorm(double x, int upper) 00036 { 00037 double ret, z, y; 00038 int up; 00039 00040 up = upper; 00041 z = x; 00042 00043 if (x < 0.0) { 00044 up = up == 0 ? 1 : 0; 00045 z = -x; 00046 } 00047 00048 if (!(z <= LTONE || (up == 1 && z <= UTZERO))) 00049 ret = 0.0; 00050 else { 00051 y = 0.5 * z * z; 00052 if (z <= 1.28) 00053 ret = 0.5 - z * (0.398942280444 - 0.399903438504 * y / 00054 (y + 5.75885480458 - 29.8213557808 / 00055 (y + 2.62433121679 + 48.6959930692 / 00056 (y + 5.92885724438)))); 00057 else 00058 ret = 0.398942280385 * exp(-y) / 00059 (z - 3.8052e-8 + 1.00000615302 / 00060 (z + 3.98064794e-4 + 1.98615381364 / 00061 (z - 0.151679116635 + 5.29330324926 / 00062 (z + 4.8385912808 - 15.1508972451 / 00063 (z + 0.742380924027 + 30.789933034 / 00064 (z + 3.99019417011)))))); 00065 } 00066 00067 if (up == 0) 00068 ret = 1.0 - ret; 00069 00070 return ret; 00071 }