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 double *watson_u2_exp(double *x, int n) 00008 { 00009 double *xcopy, mean = 0.0, zbar = 0.0, fn2, fx, sum4 = 0.0; 00010 static double y[2]; 00011 int i; 00012 00013 if ((xcopy = (double *)malloc(n * sizeof(double))) == NULL) { 00014 fprintf(stderr, "Memory error in watson_u2_exp\n"); 00015 exit(EXIT_FAILURE); 00016 } 00017 00018 for (i = 0; i < n; ++i) { 00019 xcopy[i] = x[i]; 00020 mean += x[i]; 00021 } 00022 mean /= n; 00023 00024 qsort(xcopy, n, sizeof(double), dcmp); 00025 00026 for (i = 0; i < n; ++i) { 00027 fx = 1 - exp(-xcopy[i] / mean); 00028 if (fx <= 1e-5) 00029 fx = 1e-5; 00030 00031 if (fx >= 0.99999) 00032 fx = 0.99999; 00033 00034 /* sum3 += 2 * (i + 1) * (log (fx) + log (1.0 - fx[n - i - 1])); */ 00035 fn2 = (2.0 * i + 1.0) / (2.0 * n); 00036 sum4 += (fx - fn2) * (fx - fn2); 00037 fn2 = (2.0 * (i + 1) - 1.0) / (2.0 * n); 00038 zbar += fx; 00039 } 00040 00041 zbar /= n; 00042 y[0] = (1.0 / (n * 12) + sum4) - n * (zbar - .5) * (zbar - .5); 00043 y[0] *= 1.0 + 0.16 / n; 00044 00045 #ifdef NOISY 00046 fprintf(stdout, " TEST19 WU2(E) =%10.4f\n", y[0]); 00047 #endif /* NOISY */ 00048 00049 free(xcopy); 00050 00051 return y; 00052 }