GRASS Programmer's Manual
6.4.2(2012)
|
00001 #include <math.h> 00002 00003 00004 /*- 00005 * SUBROUTINE NORMP(Z, P, Q, PDF) 00006 * 00007 * Normal distribution probabilities accurate to 1.e-15. 00008 * Z = no. of standard deviations from the mean. 00009 * P, Q = probabilities to the left & right of Z. P + Q = 1. 00010 * PDF = the probability density. 00011 * 00012 * Based upon algorithm 5666 for the error function, from: 00013 * Hart, J.F. et al, 'Computer Approximations', Wiley 1968 00014 * 00015 * Programmer: Alan Miller 00016 * 00017 * Latest revision - 30 March 1986 00018 * 00019 */ 00020 00021 /* Conversion to C by James Darrell McCauley, 24 September 1994 */ 00022 00023 double normp(double z) 00024 { 00025 static double p[7] = { 220.2068679123761, 221.2135961699311, 00026 112.079291497870, 33.91286607838300, 6.37396220353165, 00027 0.7003830644436881, 0.352624965998910e-1 00028 }; 00029 static double q[8] = { 440.4137358247522, 793.8265125199484, 00030 637.3336333788311, 296.5642487796737, 86.78073220294608, 00031 16.06417757920695, 1.755667163182642, 0.8838834764831844e-1 00032 }; 00033 static double cutoff = 7.071, root2pi = 2.506628274631001; 00034 double zabs, expntl; 00035 double pee, queue, pdf; 00036 00037 zabs = fabs(z); 00038 00039 if (zabs > 37.0) { 00040 pdf = 0.0; 00041 if (z > 0.0) { 00042 pee = 1.0; 00043 queue = 0.0; 00044 } 00045 else { 00046 pee = 0.0; 00047 queue = 1.0; 00048 } 00049 00050 return pee; 00051 } 00052 00053 expntl = exp(-0.5 * zabs * zabs); 00054 pdf = expntl / root2pi; 00055 00056 if (zabs < cutoff) 00057 pee = expntl * ((((((p[6] * zabs + p[5]) * zabs + p[4]) 00058 * zabs + p[3]) * zabs + p[2]) * zabs + 00059 p[1]) * zabs + p[0]) 00060 / (((((((q[7] * zabs + q[6]) * zabs + q[5]) * zabs + q[4]) 00061 * zabs + q[3]) * zabs + q[2]) * zabs + q[1]) * zabs + q[0]); 00062 else 00063 pee = pdf / (zabs + 1.0 / (zabs + 2.0 / (zabs + 3.0 / (zabs + 4.0 / 00064 (zabs + 00065 0.65))))); 00066 00067 if (z < 0.0) 00068 queue = 1.0 - pee; 00069 else { 00070 queue = pee; 00071 pee = 1.0 - queue; 00072 } 00073 00074 return pee; 00075 }