GRASS Programmer's Manual
6.4.2(2012)
|
00001 #include <math.h> 00002 #include <stdio.h> 00003 #include <grass/cdhc.h> 00004 #include "local_proto.h" 00005 00006 00007 /* Local function prototypes */ 00008 static double poly(double c[], int nord, double x); 00009 00010 00011 /*-Algorithm AS 181 00012 * by J.P. Royston, 1982. 00013 * Applied Statistics 31(2):176-180 00014 * 00015 * Translation to C by James Darrell McCauley, mccauley@ecn.purdue.edu. 00016 * 00017 * Calculates Shapiro and Wilk's W statistic and its sig. level 00018 * 00019 * Originally used: 00020 * Auxiliary routines required: ALNORM = algorithm AS 66 and NSCOR2 00021 * from AS 177. 00022 00023 * Note: ppnd() from as66 was replaced with ppnd16() from as241. 00024 */ 00025 void wext(double x[], int n, double ssq, double a[], int n2, double eps, 00026 double *w, double *pw, int *ifault) 00027 { 00028 double eu3, lamda, ybar, sdy, al, un, ww, y, z; 00029 int i, j, n3, nc; 00030 static double wa[3] = { 0.118898, 0.133414, 0.327907 }; 00031 static double wb[4] = { -0.37542, -0.492145, -1.124332, -0.199422 }; 00032 static double wc[4] = { -3.15805, 0.729399, 3.01855, 1.558776 }; 00033 static double wd[6] = { 0.480385, 0.318828, 0.0, -0.0241665, 0.00879701, 00034 0.002989646 00035 }; 00036 static double we[6] = { -1.91487, -1.37888, -0.04183209, 0.1066339, 00037 -0.03513666, -0.01504614 00038 }; 00039 static double wf[7] = { -3.73538, -1.015807, -0.331885, 0.1773538, 00040 -0.01638782, -0.03215018, 0.003852646 00041 }; 00042 static double unl[3] = { -3.8, -3.0, -1.0 }; 00043 static double unh[3] = { 8.6, 5.8, 5.4 }; 00044 static int nc1[3] = { 5, 5, 5 }; 00045 static int nc2[3] = { 3, 4, 5 }; 00046 double c[5]; 00047 int upper = 1; 00048 static double pi6 = 1.90985932, stqr = 1.04719755; 00049 static double zero = 0.0, tqr = 0.75, one = 1.0; 00050 static double onept4 = 1.4, three = 3.0, five = 5.0; 00051 static double c1[5][3] = { 00052 {-1.26233, -2.28135, -3.30623}, 00053 {1.87969, 2.26186, 2.76287}, 00054 {0.0649583, 0.0, -0.83484}, 00055 {-0.0475604, 0.0, 1.20857}, 00056 {-0.0139682, -0.00865763, -0.507590} 00057 }; 00058 static double c2[5][3] = { 00059 {-0.287696, -1.63638, -5.991908}, 00060 {1.78953, 5.60924, 21.04575}, 00061 {-0.180114, -3.63738, -24.58061}, 00062 {0.0, 1.08439, 13.78661}, 00063 {0.0, 0.0, -2.835295} 00064 }; 00065 00066 *ifault = 1; 00067 00068 *pw = one; 00069 *w = one; 00070 00071 if (n <= 2) 00072 return; 00073 00074 *ifault = 3; 00075 if (n / 2 != n2) 00076 return; 00077 00078 *ifault = 2; 00079 if (n > 2000) 00080 return; 00081 00082 *ifault = 0; 00083 i = n - 1; 00084 00085 for (*w = 0.0, j = 0; j < n2; ++j) 00086 *w += a[j] * (x[i--] - x[j]); 00087 00088 *w *= *w / ssq; 00089 if (*w > one) { 00090 *w = one; 00091 00092 return; 00093 } 00094 else if (n > 6) { /* Get significance level of W */ 00095 /* 00096 * N between 7 and 2000 ... Transform W to Y, get mean and sd, 00097 * standardize and get significance level 00098 */ 00099 00100 if (n <= 20) { 00101 al = log((double)n) - three; 00102 lamda = poly(wa, 3, al); 00103 ybar = exp(poly(wb, 4, al)); 00104 sdy = exp(poly(wc, 4, al)); 00105 } 00106 else { 00107 al = log((double)n) - five; 00108 lamda = poly(wd, 6, al); 00109 ybar = exp(poly(we, 6, al)); 00110 sdy = exp(poly(wf, 7, al)); 00111 } 00112 00113 y = pow(one - *w, lamda); 00114 z = (y - ybar) / sdy; 00115 *pw = alnorm(z, upper); 00116 00117 return; 00118 } 00119 else { 00120 /* Deal with N less than 7 (Exact significance level for N = 3). */ 00121 if (*w >= eps) { 00122 ww = *w; 00123 if (*w >= eps) { 00124 ww = *w; 00125 if (n == 3) { 00126 *pw = pi6 * (atan(sqrt(ww / (one - ww))) - stqr); 00127 00128 return; 00129 } 00130 00131 un = log((*w - eps) / (one - *w)); 00132 n3 = n - 3; 00133 if (un >= unl[n3 - 1]) { 00134 if (un <= onept4) { 00135 nc = nc1[n3 - 1]; 00136 00137 for (i = 0; i < nc; ++i) 00138 c[i] = c1[i][n3 - 1]; 00139 00140 eu3 = exp(poly(c, nc, un)); 00141 } 00142 else { 00143 if (un > unh[n3 - 1]) 00144 return; 00145 00146 nc = nc2[n3 - 1]; 00147 00148 for (i = 0; i < nc; ++i) 00149 c[i] = c2[i][n3 - 1]; 00150 00151 un = log(un); /*alog */ 00152 eu3 = exp(exp(poly(c, nc, un))); 00153 } 00154 ww = (eu3 + tqr) / (one + eu3); 00155 *pw = pi6 * (atan(sqrt(ww / (one - ww))) - stqr); 00156 00157 return; 00158 } 00159 } 00160 } 00161 *pw = zero; 00162 00163 return; 00164 } 00165 00166 return; 00167 } 00168 00169 00170 /* 00171 * Algorithm AS 181.1 Appl. Statist. (1982) Vol. 31, No. 2 00172 * 00173 * Obtain array A of weights for calculating W 00174 */ 00175 void wcoef(double a[], int n, int n2, double *eps, int *ifault) 00176 { 00177 static double c4[2] = { 0.6869, 0.1678 }; 00178 static double c5[2] = { 0.6647, 0.2412 }; 00179 static double c6[3] = { 0.6431, 0.2806, 0.0875 }; 00180 static double rsqrt2 = 0.70710678; 00181 double a1star, a1sq, sastar, an; 00182 int j; 00183 00184 *ifault = 1; 00185 if (n <= 2) 00186 return; 00187 00188 *ifault = 3; 00189 if (n / 2 != n2) 00190 return; 00191 00192 *ifault = 2; 00193 if (n > 2000) 00194 return; 00195 00196 *ifault = 0; 00197 if (n > 6) { 00198 /* Calculate rankits using approximate function nscor2(). (AS177) */ 00199 nscor2(a, n, n2, ifault); 00200 00201 for (sastar = 0.0, j = 1; j < n2; ++j) 00202 sastar += a[j] * a[j]; 00203 00204 sastar *= 8.0; 00205 00206 an = n; 00207 if (n <= 20) 00208 an--; 00209 a1sq = exp(log(6.0 * an + 7.0) - log(6.0 * an + 13.0) 00210 + 0.5 * (1.0 + (an - 2.0) * log(an + 1.0) - (an - 1.0) 00211 * log(an + 2.0))); 00212 a1star = sastar / (1.0 / a1sq - 2.0); 00213 sastar = sqrt(sastar + 2.0 * a1star); 00214 a[0] = sqrt(a1star) / sastar; 00215 00216 for (j = 1; j < n2; ++j) 00217 a[j] = 2.0 * a[j] / sastar; 00218 } 00219 else { 00220 /* Use exact values for weights */ 00221 00222 a[0] = rsqrt2; 00223 if (n != 3) { 00224 if (n - 3 == 3) 00225 for (j = 0; j < 3; ++j) 00226 a[j] = c6[j]; 00227 else if (n - 3 == 2) 00228 for (j = 0; j < 2; ++j) 00229 a[j] = c5[j]; 00230 else 00231 for (j = 0; j < 2; ++j) 00232 a[j] = c4[j]; 00233 00234 /*- 00235 goto (40,50,60), n3 00236 40 do 45 j = 1,2 00237 45 a(j) = c4(j) 00238 goto 70 00239 50 do 55 j = 1,2 00240 55 a(j) = c5(j) 00241 goto 70 00242 60 do 65 j = 1,3 00243 65 a(j) = c6(j) 00244 */ 00245 } 00246 } 00247 00248 /* Calculate the minimum possible value of W */ 00249 *eps = a[0] * a[0] / (1.0 - 1.0 / (double)n); 00250 00251 return; 00252 } 00253 00254 00255 /* 00256 * Algorithm AS 181.2 Appl. Statist. (1982) Vol. 31, No. 2 00257 * 00258 * Calculates the algebraic polynomial of order nored-1 with array of 00259 * coefficients c. Zero order coefficient is c(1) 00260 */ 00261 static double poly(double c[], int nord, double x) 00262 { 00263 double p; 00264 int n2, i, j; 00265 00266 if (nord == 1) 00267 return c[0]; 00268 00269 p = x * c[nord - 1]; 00270 00271 if (nord != 2) { 00272 n2 = nord - 2; 00273 j = n2; 00274 00275 for (i = 0; i < n2; ++i) 00276 p = (p + c[j--]) * x; 00277 } 00278 00279 return c[0] + p; 00280 } 00281 00282 00283 /* 00284 * AS R63 Appl. Statist. (1986) Vol. 35, No.2 00285 * 00286 * A remark on AS 181 00287 * 00288 * Calculates Sheppard corrected version of W test. 00289 * 00290 * Auxiliary functions required: ALNORM = algorithm AS 66, and PPND = 00291 * algorithm AS 111 (or PPND7 from AS 241). 00292 */ 00293 void wgp(double x[], int n, double ssq, double gp, double h, double a[], 00294 int n2, double eps, double w, double u, double p, int *ifault) 00295 { 00296 double zbar, zsd, an1, hh; 00297 00298 zbar = 0.0; 00299 zsd = 1.0; 00300 *ifault = 1; 00301 00302 if (n < 7) 00303 return; 00304 00305 if (gp > 0.0) { /* No correction applied if gp=0. */ 00306 an1 = (double)(n - 1); 00307 /* correct ssq and find standardized grouping interval (h) */ 00308 ssq = ssq - an1 * gp * gp / 12.0; 00309 h = gp / sqrt(ssq / an1); 00310 *ifault = 4; 00311 00312 if (h > 1.5) 00313 return; 00314 } 00315 wext(x, n, ssq, a, n2, eps, &w, &p, ifault); 00316 00317 if (*ifault != 0) 00318 return; 00319 00320 if (!(p > 0.0 && p < 1.0)) { 00321 u = 5.0 - 10.0 * p; 00322 00323 return; 00324 } 00325 00326 if (gp > 0.0) { 00327 /* correct u for grouping interval (n<=100 and n>100 separately) */ 00328 hh = sqrt(h); 00329 if (n <= 100) { 00330 zbar = -h * (1.07457 + hh * (-2.8185 + hh * 1.8898)); 00331 zsd = 1.0 + h * (0.50933 + hh * (-0.98305 + hh * 0.7408)); 00332 } 00333 else { 00334 zbar = -h * (0.96436 + hh * (-2.1300 + hh * 1.3196)); 00335 zsd = 1.0 + h * (0.2579 + h * 0.15225); 00336 } 00337 } 00338 00339 /* ppnd is AS 111 (Beasley and Springer, 1977) */ 00340 u = (-ppnd16(p) - zbar) / zsd; 00341 00342 /* alnorm is AS 66 (Hill, 1973) */ 00343 p = alnorm(u, 1); 00344 00345 return; 00346 }