GRASS Programmer's Manual  6.4.2(2012)
as181.c
Go to the documentation of this file.
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 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines