GRASS Programmer's Manual  6.4.2(2012)
durbins.c
Go to the documentation of this file.
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 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines