GRASS GIS 8 Programmer's Manual 8.3.2(2024)-exported
Loading...
Searching...
No Matches
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 */
7static 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 */
23void 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 */
166void 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 */
252static 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 */
283void 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