GRASS GIS 8 Programmer's Manual 8.2.1RC1(2022)-exported
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 */
8static 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 */
25void 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 */
175void 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 */
261static 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 */
293void 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}
void Cdhc_nscor2(double s[], int n, int n2, int *ifault)
Definition: as177.c:111
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 wcoef(double a[], int n, int n2, double *eps, int *ifault)
Definition: as181.c:175
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
double ppnd16(double p)
Definition: as241.c:90
double Cdhc_alnorm(double x, int upper)
Definition: as66.c:35
#define x