-
-
Notifications
You must be signed in to change notification settings - Fork 288
/
royston.c
65 lines (52 loc) · 1.27 KB
/
royston.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "local_proto.h"
/*-
* driver program for AS 181: Royston's extension of the Shapiro-Wilk
* W statistic to n=2000
* needs as181.c as177.c as241.c Cdhc_dcmp.c as66.c
*/
double *Cdhc_royston(double *x, int n)
{
static double y[2];
double *a, eps, w, pw, mean = 0, ssq = 0, *xcopy;
int i, ifault, n2;
n2 = (int)floor((double)n / 2);
#ifndef lint
if ((a = (double *)malloc(n2 * sizeof(double))) == NULL) {
fprintf(stderr, "Memory error in royston\n");
exit(EXIT_FAILURE);
}
if ((xcopy = (double *)malloc(n * sizeof(double))) == NULL) {
fprintf(stderr, "Memory error in royston\n");
exit(EXIT_FAILURE);
}
#endif /* lint */
for (i = 0; i < n; ++i) {
xcopy[i] = x[i];
mean += x[i];
}
mean /= n;
qsort(xcopy, n, sizeof(double), Cdhc_dcmp);
for (i = 0; i < n; ++i)
ssq += (mean - x[i]) * (mean - x[i]);
wcoef(a, n, n2, &eps, &ifault);
if (ifault == 0)
wext(xcopy, n, ssq, a, n2, eps, &w, &pw, &ifault);
else {
fprintf(stderr, "Error in wcoef()\n");
return (double *)NULL;
}
if (ifault == 0) {
y[0] = w;
y[1] = pw;
}
else {
fprintf(stderr, "Error in wcoef()\n");
return (double *)NULL;
}
free(a);
free(xcopy);
return y;
}