GRASS GIS 8 Programmer's Manual 8.2.1RC1(2022)-exported
durbins.c
Go to the documentation of this file.
1#include <stdio.h>
2#include <stdlib.h>
3#include <math.h>
4#include "local_proto.h"
5
6
7/* could probably use some cleanup/optimization */
8double *Cdhc_durbins_exact(double *x, int n)
9{
10 static double y[2];
11 double *xcopy, sumx = 0.0, sumx2 = 0.0, s2, *b, *c, *g, *z, sqrt2;
12 int i, j;
13
14 if ((b = (double *)malloc(n * sizeof(double))) == NULL) {
15 fprintf(stderr, "Memory error in Cdhc_durbins_exact\n");
16 exit(EXIT_FAILURE);
17 }
18 if ((c = (double *)malloc((n + 1) * sizeof(double))) == NULL) {
19 fprintf(stderr, "Memory error in Cdhc_durbins_exact\n");
20 exit(EXIT_FAILURE);
21 }
22 if ((g = (double *)malloc((n + 1) * sizeof(double))) == NULL) {
23 fprintf(stderr, "Memory error in Cdhc_durbins_exact\n");
24 exit(EXIT_FAILURE);
25 }
26 if ((z = (double *)malloc(n * sizeof(double))) == NULL) {
27 fprintf(stderr, "Memory error in Cdhc_durbins_exact\n");
28 exit(EXIT_FAILURE);
29 }
30 if ((xcopy = (double *)malloc(n * sizeof(double))) == NULL) {
31 fprintf(stderr, "Memory error in Cdhc_durbins_exact\n");
32 exit(EXIT_FAILURE);
33 }
34
35 sqrt2 = sqrt((double)2.0);
36 for (i = 0; i < n; ++i) {
37 xcopy[i] = x[i];
38 sumx += x[i];
39 sumx2 += x[i] * x[i];
40 }
41
42 s2 = sqrt((sumx2 - sumx * sumx / n) / (n - 1));
43 for (i = 0; i < n; ++i) {
44 xcopy[i] = (xcopy[i] - sumx / n) / s2;
45 b[i] = 0.5 + Cdhc_normp(xcopy[i] / sqrt2) / 2.0;
46 }
47
48 qsort(b, n, sizeof(double), Cdhc_dcmp);
49
50 for (i = 1; i < n; ++i)
51 c[i] = b[i] - b[i - 1];
52
53 c[0] = b[0];
54 c[n] = 1 - b[n - 1];
55
56 qsort(c, n + 1, sizeof(double), Cdhc_dcmp);
57
58 for (j = 1; j <= n; ++j)
59 g[j] = (n + 1 - j) * (c[j] - c[j - 1]);
60
61 g[0] = (n + 1) * c[0];
62 g[n] = c[n] - c[n - 1];
63
64 for (i = 0; i < n; ++i) {
65 z[i] = 0.0;
66 for (j = 0; j <= i; ++j)
67 z[i] += g[j];
68
69 z[i] = (i + 1.0) / n - z[i];
70 }
71
72 qsort(z, n, sizeof(double), Cdhc_dcmp);
73
74 y[0] = z[n - 1];
75 y[1] = sqrt((double)n) * z[n - 1];
76
77#ifdef NOISY
78 fprintf(stdout, " TEST7 DRB(N) =%10.4f\n", y[0]);
79#endif /* NOISY */
80
81 free(b);
82 free(c);
83 free(g);
84 free(xcopy);
85 free(z);
86
87 return y;
88}
#define NULL
Definition: ccmath.h:32
int Cdhc_dcmp(const void *i, const void *j)
Definition: dcmp.c:1
double b
double * Cdhc_durbins_exact(double *x, int n)
Definition: durbins.c:8
float g
Definition: named_colr.c:8
double Cdhc_normp(double z)
Definition: normp.c:23
#define x