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