GRASS Programmer's Manual
6.4.2(2012)
|
00001 #include <stdio.h> 00002 #include <stdlib.h> 00003 #include <math.h> 00004 #include "local_proto.h" 00005 00006 00007 /* could probably use some cleanup/optimization */ 00008 double *durbins_exact(double *x, int n) 00009 { 00010 static double y[2]; 00011 double *xcopy, sumx = 0.0, sumx2 = 0.0, s2, *b, *c, *g, *z, sqrt2; 00012 int i, j; 00013 00014 if ((b = (double *)malloc(n * sizeof(double))) == NULL) { 00015 fprintf(stderr, "Memory error in durbins_exact\n"); 00016 exit(EXIT_FAILURE); 00017 } 00018 if ((c = (double *)malloc((n + 1) * sizeof(double))) == NULL) { 00019 fprintf(stderr, "Memory error in durbins_exact\n"); 00020 exit(EXIT_FAILURE); 00021 } 00022 if ((g = (double *)malloc((n + 1) * sizeof(double))) == NULL) { 00023 fprintf(stderr, "Memory error in durbins_exact\n"); 00024 exit(EXIT_FAILURE); 00025 } 00026 if ((z = (double *)malloc(n * sizeof(double))) == NULL) { 00027 fprintf(stderr, "Memory error in durbins_exact\n"); 00028 exit(EXIT_FAILURE); 00029 } 00030 if ((xcopy = (double *)malloc(n * sizeof(double))) == NULL) { 00031 fprintf(stderr, "Memory error in durbins_exact\n"); 00032 exit(EXIT_FAILURE); 00033 } 00034 00035 sqrt2 = sqrt((double)2.0); 00036 for (i = 0; i < n; ++i) { 00037 xcopy[i] = x[i]; 00038 sumx += x[i]; 00039 sumx2 += x[i] * x[i]; 00040 } 00041 00042 s2 = sqrt((sumx2 - sumx * sumx / n) / (n - 1)); 00043 for (i = 0; i < n; ++i) { 00044 xcopy[i] = (xcopy[i] - sumx / n) / s2; 00045 b[i] = 0.5 + normp(xcopy[i] / sqrt2) / 2.0; 00046 } 00047 00048 qsort(b, n, sizeof(double), dcmp); 00049 00050 for (i = 1; i < n; ++i) 00051 c[i] = b[i] - b[i - 1]; 00052 00053 c[0] = b[0]; 00054 c[n] = 1 - b[n - 1]; 00055 00056 qsort(c, n + 1, sizeof(double), dcmp); 00057 00058 for (j = 1; j <= n; ++j) 00059 g[j] = (n + 1 - j) * (c[j] - c[j - 1]); 00060 00061 g[0] = (n + 1) * c[0]; 00062 g[n] = c[n] - c[n - 1]; 00063 00064 for (i = 0; i < n; ++i) { 00065 z[i] = 0.0; 00066 for (j = 0; j <= i; ++j) 00067 z[i] += g[j]; 00068 00069 z[i] = (i + 1.0) / n - z[i]; 00070 } 00071 00072 qsort(z, n, sizeof(double), dcmp); 00073 00074 y[0] = z[n - 1]; 00075 y[1] = sqrt((double)n) * z[n - 1]; 00076 00077 #ifdef NOISY 00078 fprintf(stdout, " TEST7 DRB(N) =%10.4f\n", y[0]); 00079 #endif /* NOISY */ 00080 00081 free(b); 00082 free(c); 00083 free(g); 00084 free(xcopy); 00085 free(z); 00086 00087 return y; 00088 }