GRASS Programmer's Manual  6.5.svn(2014)-r66266
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
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 */
8 double *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 durbins_exact\n");
16  exit(EXIT_FAILURE);
17  }
18  if ((c = (double *)malloc((n + 1) * sizeof(double))) == NULL) {
19  fprintf(stderr, "Memory error in durbins_exact\n");
20  exit(EXIT_FAILURE);
21  }
22  if ((g = (double *)malloc((n + 1) * sizeof(double))) == NULL) {
23  fprintf(stderr, "Memory error in durbins_exact\n");
24  exit(EXIT_FAILURE);
25  }
26  if ((z = (double *)malloc(n * sizeof(double))) == NULL) {
27  fprintf(stderr, "Memory error in durbins_exact\n");
28  exit(EXIT_FAILURE);
29  }
30  if ((xcopy = (double *)malloc(n * sizeof(double))) == NULL) {
31  fprintf(stderr, "Memory error in 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 + normp(xcopy[i] / sqrt2) / 2.0;
46  }
47 
48  qsort(b, n, sizeof(double), 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), 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), 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 }
int dcmp(const void *i, const void *j)
Definition: dcmp.c:1
float b
Definition: named_colr.c:8
double * durbins_exact(double *x, int n)
Definition: durbins.c:8
int y
Definition: plot.c:34
void * malloc(YYSIZE_T)
float g
Definition: named_colr.c:8
double normp(double)
Definition: normp.c:23
return NULL
Definition: dbfopen.c:1394
void free(void *)
int n
Definition: dataquad.c:291