Changeset 36290

Timestamp:

Nov 20, 2013, 3:21:02 PM (13 years ago)

Author:

watersc1

Message:

This should resolve the issues with the fitted statistics. Briefly,
the prior change to enforce 5-sigma clipping thresholds altered the
effect of the choice of binsize. Because of this, we we making the
histogram for the fitting data with binsizes far too small to be
appropriate. This led to very noisy histograms, and poor centering
for the fit, and so very bad fits, with very bad statistics.

I've changed the binsize calculation to attempt to keep the 2-sigma
bin filled with ~50 data samples. This seems to produce a usable
histogram, and the fitted statistics match better. There is a slight
reduction in the mean bias from the random data test, and some
improvement in the associated sigma values. This appears to be the
same in the ppSim test.

File:

• trunk/psLib/src/math/psStats.c (modified) (10 diffs)

trunk/psLib/src/math/psStats.c

@@ -1049 +1049 @@
         
 #if (CZW)
-        printf("CZW: bad sigma?: %f %f  %f %f  %f %f  %f %f %f  %f\n",
-               binH2F32,binL2F32,binHiF32,binLoF32,binH4F32,binL4F32,
-               sigma1,sigma2,sigma4,sigma);
+        //      printf("CZW: bad sigma?: %f %f  %f %f  %f %f  %f %f %f  %f\n",
+        //             binH2F32,binL2F32,binHiF32,binLoF32,binH4F32,binL4F32,
+        //             sigma1,sigma2,sigma4,sigma);
         
-        printf("CZW (%d): median %f sigma %f delta: %f \n\t %f %f %f %f %f %f %f \n\t %f %f %f %f %f %f %f\n",
+        printf("CZW Robust (%d): median %f sigma %f delta: %f \n\t %f %f %f %f %f %f %f \n\t %f %f %f %f %f %f %f\n",
                iterate,
-           stats->robustMedian,stats->robustStdev,
-           fabs(cumulative->bounds->data.F32[binMedian] - cumulative->bounds->data.F32[binMedian + 1]),
-
-           cumulative->bounds->data.F32[binMedian-3],cumulative->bounds->data.F32[binMedian-2],
-           cumulative->bounds->data.F32[binMedian-1],
-           cumulative->bounds->data.F32[binMedian],
-           cumulative->bounds->data.F32[binMedian+1],
-           cumulative->bounds->data.F32[binMedian+2],cumulative->bounds->data.F32[binMedian+3],
-
-           cumulative->nums->data.F32[binMedian-3],cumulative->nums->data.F32[binMedian-2],
-           cumulative->nums->data.F32[binMedian-1],
-           cumulative->nums->data.F32[binMedian],
-           cumulative->nums->data.F32[binMedian+1],
-           cumulative->nums->data.F32[binMedian+2],cumulative->nums->data.F32[binMedian+3]);
+               stats->robustMedian,stats->robustStdev,
+               fabs(cumulative->bounds->data.F32[binMedian] - cumulative->bounds->data.F32[binMedian + 1]),
+               
+               cumulative->bounds->data.F32[binMedian-3],cumulative->bounds->data.F32[binMedian-2],
+               cumulative->bounds->data.F32[binMedian-1],
+               cumulative->bounds->data.F32[binMedian],
+               cumulative->bounds->data.F32[binMedian+1],
+               cumulative->bounds->data.F32[binMedian+2],cumulative->bounds->data.F32[binMedian+3],
+               
+               cumulative->nums->data.F32[binMedian-3],cumulative->nums->data.F32[binMedian-2],
+               cumulative->nums->data.F32[binMedian-1],
+               cumulative->nums->data.F32[binMedian],
+               cumulative->nums->data.F32[binMedian+1],
+               cumulative->nums->data.F32[binMedian+2],cumulative->nums->data.F32[binMedian+3]);
         //      PS_VECTOR_PRINT_F32(histogram->bounds);
         //      PS_VECTOR_PRINT_F32(histogram->nums);
@@ -1262 +1262 @@
             // set roughly so that the lowest bins have about 2 cnts
             // Nsmallest ~ N50 / (4*dN))
-            psF32 dN = PS_MAX (1, PS_MIN (4, stats->robustN50 / 8)); 
-           binSize = guessStdev / dN;
+          //            psF32 dN = PS_MAX (1, PS_MIN (4, stats->robustN50 / 8));
+
+          // CZW 2013-11-20: We know that the histogram is going to be basically Gaussian.
+          // Furthermore, we only use the inner +/- 2 sigma parts.  Therefore, define the
+          // binsize such that the bin at 2 sigma contains ~50 points (S/N ~ 7).  robustN50
+          // contains half the total points, so 2 * robustN50 / 50 is the fraction of all
+          // points in the 2 sigma bin.  Dance the erf() relations around, and it looks like
+          // there's a factor of about 1/20 to include.  Keep the PS_MAX to ensure we never bin
+          // wider than 1 sigma when the number of points is small.
+          psF32 dN = PS_MAX(1, (stats->robustN50 / 500.0));
+          binSize = guessStdev / dN;
         }
 
@@ -1291 +1300 @@
         // We no longer want to reset the binSize here, as it can cause odd things.  Better to select
         // a number of bins, and then set the min/max values to put those bins sanely around the mean.
-        long numBins = PS_MAX (50, PS_MIN (10000, (max - min) / binSize));
+        //        long numBins = PS_MAX (50, PS_MIN (10000, (max - min) / binSize));
         //        binSize = (max - min) / (float) numBins;
         psTrace(TRACE, 6, "The new min/max values are (%f, %f).\n", min, max);
         psTrace(TRACE, 6, "The new bin size is %f.\n", binSize);
-        psTrace(TRACE, 6, "The numBins is %ld\n", numBins);
+        //        psTrace(TRACE, 6, "The numBins is %ld\n", numBins);
 
 
@@ -1305 +1314 @@
           max = guessMean + FITTED_CLIPPING_NUM * guessStdev;
         }
-        
+        long numBins = PS_MAX (50, PS_MIN (10000, (max - min) / binSize));
+        if (CZW) { printf("I've clipped: %f %f => %f %f ; %f %ld\n",guessMean,guessStdev,min,max,binSize,stats->robustN50); }
         psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers)
         if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) {
@@ -1357 +1367 @@
             }
             psTrace (TRACE, 6, "(%f = %.0f) ", histogram->bounds->data.F32[i], histogram->nums->data.F32[i]);
+            if (CZW) { printf("CENTERED_HISTOGRAM: %f %f\n",
+                              PS_BIN_MIDPOINT(histogram,i),
+                              histogram->nums->data.F32[i]); }
         }
         psTrace (TRACE, 6, "\n");
 
+        if (CZW) { printf("Bin selection done: %ld %f %f %ld %f %f %ld %f %f\n",
+                          binMin,PS_BIN_MIDPOINT(histogram,binMin),histogram->nums->data.F32[binMin],
+                          binMax,PS_BIN_MIDPOINT(histogram,binMax),histogram->nums->data.F32[binMax],
+                          binPeak,PS_BIN_MIDPOINT(histogram,binPeak),histogram->nums->data.F32[binPeak]);
+        }
+        
         // assume a reasonably well-defined gaussian-like population; run from peak out until val < 0.25*peak
         psTrace(TRACE, 6, "The clipped numBins is %ld\n", binMax - binMin);
@@ -1411 +1430 @@
             psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
             bool status = psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x);
+#if (CZW && 1)
+            printf("CZW: LowfitPoly: %f %f %f\n",poly->coeff[0],poly->coeff[1],poly->coeff[2]);
+            for (long i = 0; i < x->n; i++) {
+              printf("CZW: Lowfit: %d %ld %f %f %f\n",
+                     status,i,x->data.F32[i],y->data.F32[i],
+                     poly->coeff[0] + poly->coeff[1] * x->data.F32[i] +
+                     poly->coeff[2] * pow(x->data.F32[i],2));
+            }
+#endif
             psFree(x);
             psFree(y);
@@ -1417 +1445 @@
                 psErrorClear();
                 COUNT_WARNING(10, 100, "Failed to fit a gaussian to the robust histogram.\n");
-
                 psFree(poly);
                 psFree(histogram);
@@ -1426 +1453 @@
             if (poly->coeff[2] >= 0.0) {
                 COUNT_WARNING(10, 100, "Failed parabolic fit: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]);
+
                 psFree(poly);
                 psFree(histogram);
@@ -1478 +1506 @@
                 }
             }
+
             psTrace(TRACE, 6, "Lower bound for symmetric range: %f (%ld)\n",
                     PS_BIN_MIDPOINT(histogram, binS), binS);
@@ -1499 +1528 @@
             psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
             bool status = psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x);
-#if (CZW && 0)
+#if (CZW && 1)
+            printf("CZW: FullfitPoly: %f %f %f\n",poly->coeff[0],poly->coeff[1],poly->coeff[2]);
             for (long i = 0; i < x->n; i++) {
-              printf("CZW: Dcheck: %ld %f %f %f\n",
-                     i,x->data.F32[i],y->data.F32[i],
+              printf("CZW: Fullfit: %d %ld %f %f %f\n",
+                     status,i,x->data.F32[i],y->data.F32[i],
                      poly->coeff[0] + poly->coeff[1] * x->data.F32[i] +
                      poly->coeff[2] * pow(x->data.F32[i],2));