GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-6c790bf5c0
as181.c
Go to the documentation of this file.
1 #include <math.h>
2 #include <stdio.h>
3 #include <grass/cdhc.h>
4 #include "local_proto.h"
5 
6 /* Local function prototypes */
7 static double poly(double c[], int nord, double x);
8 
9 /*-Algorithm AS 181
10  * by J.P. Royston, 1982.
11  * Applied Statistics 31(2):176-180
12  *
13  * Translation to C by James Darrell McCauley, mccauley@ecn.purdue.edu.
14  *
15  * Calculates Shapiro and Wilk's W statistic and its sig. level
16  *
17  * Originally used:
18  * Auxiliary routines required: Cdhc_alnorm = algorithm AS 66 and Cdhc_nscor2
19  * from AS 177.
20 
21  * Note: ppnd() from as66 was replaced with ppnd16() from as241.
22  */
23 void wext(double x[], int n, double ssq, double a[], int n2, double eps,
24  double *w, double *pw, int *ifault)
25 {
26  double eu3, lamda, ybar, sdy, al, un, ww, y, z;
27  int i, j, n3, nc;
28  static double wa[3] = {0.118898, 0.133414, 0.327907};
29  static double wb[4] = {-0.37542, -0.492145, -1.124332, -0.199422};
30  static double wc[4] = {-3.15805, 0.729399, 3.01855, 1.558776};
31  static double wd[6] = {0.480385, 0.318828, 0.0,
32  -0.0241665, 0.00879701, 0.002989646};
33  static double we[6] = {-1.91487, -1.37888, -0.04183209,
34  0.1066339, -0.03513666, -0.01504614};
35  static double wf[7] = {-3.73538, -1.015807, -0.331885, 0.1773538,
36  -0.01638782, -0.03215018, 0.003852646};
37  static double unl[3] = {-3.8, -3.0, -1.0};
38  static double unh[3] = {8.6, 5.8, 5.4};
39  static int nc1[3] = {5, 5, 5};
40  static int nc2[3] = {3, 4, 5};
41  double c[5];
42  int upper = 1;
43  static double pi6 = 1.90985932, stqr = 1.04719755;
44  static double zero = 0.0, tqr = 0.75, one = 1.0;
45  static double onept4 = 1.4, three = 3.0, five = 5.0;
46 
47  static double c1[5][3] = {{-1.26233, -2.28135, -3.30623},
48  {1.87969, 2.26186, 2.76287},
49  {0.0649583, 0.0, -0.83484},
50  {-0.0475604, 0.0, 1.20857},
51  {-0.0139682, -0.00865763, -0.507590}};
52  static double c2[5][3] = {{-0.287696, -1.63638, -5.991908},
53  {1.78953, 5.60924, 21.04575},
54  {-0.180114, -3.63738, -24.58061},
55  {0.0, 1.08439, 13.78661},
56  {0.0, 0.0, -2.835295}};
57 
58  *ifault = 1;
59 
60  *pw = one;
61  *w = one;
62 
63  if (n <= 2)
64  return;
65 
66  *ifault = 3;
67  if (n / 2 != n2)
68  return;
69 
70  *ifault = 2;
71  if (n > 2000)
72  return;
73 
74  *ifault = 0;
75  i = n - 1;
76 
77  for (*w = 0.0, j = 0; j < n2; ++j)
78  *w += a[j] * (x[i--] - x[j]);
79 
80  *w *= *w / ssq;
81  if (*w > one) {
82  *w = one;
83 
84  return;
85  }
86  else if (n > 6) { /* Get significance level of W */
87  /*
88  * N between 7 and 2000 ... Transform W to Y, get mean and sd,
89  * standardize and get significance level
90  */
91 
92  if (n <= 20) {
93  al = log((double)n) - three;
94  lamda = poly(wa, 3, al);
95  ybar = exp(poly(wb, 4, al));
96  sdy = exp(poly(wc, 4, al));
97  }
98  else {
99  al = log((double)n) - five;
100  lamda = poly(wd, 6, al);
101  ybar = exp(poly(we, 6, al));
102  sdy = exp(poly(wf, 7, al));
103  }
104 
105  y = pow(one - *w, lamda);
106  z = (y - ybar) / sdy;
107  *pw = Cdhc_alnorm(z, upper);
108 
109  return;
110  }
111  else {
112  /* Deal with N less than 7 (Exact significance level for N = 3). */
113  if (*w >= eps) {
114  ww = *w;
115  if (*w >= eps) {
116  ww = *w;
117  if (n == 3) {
118  *pw = pi6 * (atan(sqrt(ww / (one - ww))) - stqr);
119 
120  return;
121  }
122 
123  un = log((*w - eps) / (one - *w));
124  n3 = n - 3;
125  if (un >= unl[n3 - 1]) {
126  if (un <= onept4) {
127  nc = nc1[n3 - 1];
128 
129  for (i = 0; i < nc; ++i)
130  c[i] = c1[i][n3 - 1];
131 
132  eu3 = exp(poly(c, nc, un));
133  }
134  else {
135  if (un > unh[n3 - 1])
136  return;
137 
138  nc = nc2[n3 - 1];
139 
140  for (i = 0; i < nc; ++i)
141  c[i] = c2[i][n3 - 1];
142 
143  un = log(un); /*alog */
144  eu3 = exp(exp(poly(c, nc, un)));
145  }
146  ww = (eu3 + tqr) / (one + eu3);
147  *pw = pi6 * (atan(sqrt(ww / (one - ww))) - stqr);
148 
149  return;
150  }
151  }
152  }
153  *pw = zero;
154 
155  return;
156  }
157 
158  return;
159 }
160 
161 /*
162  * Algorithm AS 181.1 Appl. Statist. (1982) Vol. 31, No. 2
163  *
164  * Obtain array A of weights for calculating W
165  */
166 void wcoef(double a[], int n, int n2, double *eps, int *ifault)
167 {
168  static double c4[2] = {0.6869, 0.1678};
169  static double c5[2] = {0.6647, 0.2412};
170  static double c6[3] = {0.6431, 0.2806, 0.0875};
171  static double rsqrt2 = 0.70710678;
172  double a1star, a1sq, sastar, an;
173  int j;
174 
175  *ifault = 1;
176  if (n <= 2)
177  return;
178 
179  *ifault = 3;
180  if (n / 2 != n2)
181  return;
182 
183  *ifault = 2;
184  if (n > 2000)
185  return;
186 
187  *ifault = 0;
188  if (n > 6) {
189  /* Calculate rankits using approximate function Cdhc_nscor2(). (AS177)
190  */
191  Cdhc_nscor2(a, n, n2, ifault);
192 
193  for (sastar = 0.0, j = 1; j < n2; ++j)
194  sastar += a[j] * a[j];
195 
196  sastar *= 8.0;
197 
198  an = n;
199  if (n <= 20)
200  an--;
201  a1sq = exp(log(6.0 * an + 7.0) - log(6.0 * an + 13.0) +
202  0.5 * (1.0 + (an - 2.0) * log(an + 1.0) -
203  (an - 1.0) * log(an + 2.0)));
204  a1star = sastar / (1.0 / a1sq - 2.0);
205  sastar = sqrt(sastar + 2.0 * a1star);
206  a[0] = sqrt(a1star) / sastar;
207 
208  for (j = 1; j < n2; ++j)
209  a[j] = 2.0 * a[j] / sastar;
210  }
211  else {
212  /* Use exact values for weights */
213 
214  a[0] = rsqrt2;
215  if (n != 3) {
216  if (n - 3 == 3)
217  for (j = 0; j < 3; ++j)
218  a[j] = c6[j];
219  else if (n - 3 == 2)
220  for (j = 0; j < 2; ++j)
221  a[j] = c5[j];
222  else
223  for (j = 0; j < 2; ++j)
224  a[j] = c4[j];
225 
226  /*-
227  goto (40,50,60), n3
228  40 do 45 j = 1,2
229  45 a(j) = c4(j)
230  goto 70
231  50 do 55 j = 1,2
232  55 a(j) = c5(j)
233  goto 70
234  60 do 65 j = 1,3
235  65 a(j) = c6(j)
236  */
237  }
238  }
239 
240  /* Calculate the minimum possible value of W */
241  *eps = a[0] * a[0] / (1.0 - 1.0 / (double)n);
242 
243  return;
244 }
245 
246 /*
247  * Algorithm AS 181.2 Appl. Statist. (1982) Vol. 31, No. 2
248  *
249  * Calculates the algebraic polynomial of order nored-1 with array of
250  * coefficients c. Zero order coefficient is c(1)
251  */
252 static double poly(double c[], int nord, double x)
253 {
254  double p;
255  int n2, i, j;
256 
257  if (nord == 1)
258  return c[0];
259 
260  p = x * c[nord - 1];
261 
262  if (nord != 2) {
263  n2 = nord - 2;
264  j = n2;
265 
266  for (i = 0; i < n2; ++i)
267  p = (p + c[j--]) * x;
268  }
269 
270  return c[0] + p;
271 }
272 
273 /*
274  * AS R63 Appl. Statist. (1986) Vol. 35, No.2
275  *
276  * A remark on AS 181
277  *
278  * Calculates Sheppard Cdhc_corrected version of W test.
279  *
280  * Auxiliary functions required: Cdhc_alnorm = algorithm AS 66, and PPND =
281  * algorithm AS 111 (or Cdhc_ppnd7 from AS 241).
282  */
283 void Cdhc_wgp(double x[], int n, double ssq, double gp, double h, double a[],
284  int n2, double eps, double w, double u, double p, int *ifault)
285 {
286  double zbar, zsd, an1, hh;
287 
288  zbar = 0.0;
289  zsd = 1.0;
290  *ifault = 1;
291 
292  if (n < 7)
293  return;
294 
295  if (gp > 0.0) { /* No Cdhc_correction applied if gp=0. */
296  an1 = (double)(n - 1);
297  /* Cdhc_correct ssq and find standardized grouping interval (h) */
298  ssq = ssq - an1 * gp * gp / 12.0;
299  h = gp / sqrt(ssq / an1);
300  *ifault = 4;
301 
302  if (h > 1.5)
303  return;
304  }
305  wext(x, n, ssq, a, n2, eps, &w, &p, ifault);
306 
307  if (*ifault != 0)
308  return;
309 
310  if (!(p > 0.0 && p < 1.0)) {
311  u = 5.0 - 10.0 * p;
312 
313  return;
314  }
315 
316  if (gp > 0.0) {
317  /* Cdhc_correct u for grouping interval (n<=100 and n>100 separately) */
318  hh = sqrt(h);
319  if (n <= 100) {
320  zbar = -h * (1.07457 + hh * (-2.8185 + hh * 1.8898));
321  zsd = 1.0 + h * (0.50933 + hh * (-0.98305 + hh * 0.7408));
322  }
323  else {
324  zbar = -h * (0.96436 + hh * (-2.1300 + hh * 1.3196));
325  zsd = 1.0 + h * (0.2579 + h * 0.15225);
326  }
327  }
328 
329  /* ppnd is AS 111 (Beasley and Springer, 1977) */
330  u = (-ppnd16(p) - zbar) / zsd;
331 
332  /* Cdhc_alnorm is AS 66 (Hill, 1973) */
333  p = Cdhc_alnorm(u, 1);
334 
335  return;
336 }
void Cdhc_nscor2(double s[], int n, int n2, int *ifault)
Definition: as177.c:105
void wext(double x[], int n, double ssq, double a[], int n2, double eps, double *w, double *pw, int *ifault)
Definition: as181.c:23
void wcoef(double a[], int n, int n2, double *eps, int *ifault)
Definition: as181.c:166
void Cdhc_wgp(double x[], int n, double ssq, double gp, double h, double a[], int n2, double eps, double w, double u, double p, int *ifault)
Definition: as181.c:283
double ppnd16(double p)
Definition: as241.c:86
double Cdhc_alnorm(double x, int upper)
Definition: as66.c:32
#define x