00001 #ident "$Header: /home/cvsuser/cvsroot/JSOC/proj/libs/stats/dstats2.c,v 1.4 2013/07/05 18:35:19 phil Exp $"
00002
00003 #include <stdlib.h>
00004 #include <math.h>
00005 #include <float.h>
00006 #include <string.h>
00007
00008 #define NBINS 65536
00009 #define OK 0
00010 #define TOO_FEW_GOOD_POINTS -2
00011
00012 int dstats2(int n, double arr[], double *min, double *max, double *medn,
00013 double *mean, double *sig, double *skew, double *kurt, int *ngood)
00014 {
00015 int i, bin, nv = 0;
00016 double dmin = DBL_MAX, dmax = -DBL_MAX;
00017 double s = 0.0, s2 = 0.0, s3 = 0.0, s4 = 0.0, avg, var;
00018 static int hist[NBINS];
00019 double delta;
00020
00021 *min = *max = *medn = *mean = *sig = *skew = *kurt = __builtin_nan("");
00022
00023 for (i = 0; i < n; ++i) {
00024 double t = arr[i];
00025 int fpclass = fpclassify(t);
00026 if (fpclass == FP_NAN || fpclass == FP_INFINITE)
00027 continue;
00028 if (dmin > t)
00029 dmin = t;
00030 if (dmax < t)
00031 dmax = t;
00032 s += t;
00033 ++nv;
00034 }
00035 *ngood = nv;
00036
00037 if (nv < 2) {
00038 if (nv == 1)
00039 *min = *max = *medn = *mean = s;
00040 return TOO_FEW_GOOD_POINTS;
00041 }
00042
00043 if (dmin == dmax) {
00044 *min = *max = *medn = *mean = dmin;
00045 *sig = *skew = *kurt = 0;
00046 return OK;
00047 }
00048
00049 delta = (dmax - dmin) / NBINS;
00050 memset(hist, 0, sizeof(int)*NBINS);
00051 avg = s / nv;
00052 for (i = 0; i < n; ++i) {
00053 double tmp = arr[i];
00054 double t = tmp - avg;
00055 double tt;
00056 int fpclass = fpclassify(t);
00057 if (fpclass == FP_NAN || fpclass == FP_INFINITE)
00058 continue;
00059 tt = t * t;
00060 s2 += tt;
00061 tt *= t;
00062 s3 += tt;
00063 tt *= t;
00064 s4 += tt;
00065 bin = (tmp - dmin) / delta;
00066 if (bin >= NBINS)
00067 bin = NBINS - 1;
00068 ++hist[bin];
00069 }
00070
00071 bin = 0;
00072 n = hist[bin];
00073 while (2*n < nv)
00074 n += hist[++bin];
00075 *medn = dmin + delta * (bin + 0.5);
00076
00077 *min = dmin;
00078 *max = dmax;
00079 *mean = avg;
00080 var = s2 / (nv - 1);
00081 *sig = sqrt(var);
00082 *skew = s3 / (var * *sig * nv);
00083 *kurt = s4 / (var * var * nv) - 3.0;
00084
00085 return OK;
00086 }
