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