Changeset 15126

Timestamp:

Sep 29, 2007, 1:58:55 PM (19 years ago)

Author:

eugene

Message:

for Robust Stats:

• measure +/- 0.5, 1.0, 2.0 sigma fractions to estimate sigma
• take the minimum of the three :
  • if the distribution has very wide outliers, 2.0 sigma will be biased high
  • if the distribution is bi-modal, 0.5 sigma will be biased high

(I don't currently do anything smart if the distribution is bi-modal)

for Fittest Stats (v4):

• if the fit to the parabola is insane, don't let the mean go outside the range of binMin - binMax.

File:

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

trunk/psLib/src/math/psStats.c

@@ -13 +13 @@
  * use ->min and ->max (PS_STAT_USE_RANGE)
  *
- *  @version $Revision: 1.217 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-28 00:36:30 $
+ *  @version $Revision: 1.218 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-29 23:58:55 $
  *
  *  Copyright 2006 IfA, University of Hawaii
@@ -744 +744 @@
     long binLo, binHi;
     long binL2, binH2;
+    long binL4, binH4;
     long binMedian;
 
@@ -842 +843 @@
         PS_BIN_FOR_VALUE(binL2, cumulative->nums, totalDataPoints * 0.308538f, 0);
         PS_BIN_FOR_VALUE(binH2, cumulative->nums, totalDataPoints * 0.691462f, 0);
+        PS_BIN_FOR_VALUE(binL4, cumulative->nums, totalDataPoints * 0.022481f, 0);
+        PS_BIN_FOR_VALUE(binH4, cumulative->nums, totalDataPoints * 0.977519f, 0);
 
         psTrace(TRACE, 6, "The 15.8655%% and 84.1345%% data point bins are (%ld, %ld).\n",
@@ -864 +867 @@
         // containing the value of interest.  constrain the answer to be within the current bin.
         // (extrapolation should not be needed and will result in errors)
-        float binLoF32, binHiF32, binL2F32, binH2F32;
-        PS_BIN_INTERPOLATE (binLoF32, cumulative->nums, cumulative->bounds, binL2, totalDataPoints * 0.158655f);
-        PS_BIN_INTERPOLATE (binHiF32, cumulative->nums, cumulative->bounds, binH2, totalDataPoints * 0.841345f);
+        float binLoF32, binHiF32, binL2F32, binH2F32, binL4F32, binH4F32;
+        PS_BIN_INTERPOLATE (binLoF32, cumulative->nums, cumulative->bounds, binLo, totalDataPoints * 0.158655f);
+        PS_BIN_INTERPOLATE (binHiF32, cumulative->nums, cumulative->bounds, binHi, totalDataPoints * 0.841345f);
         PS_BIN_INTERPOLATE (binL2F32, cumulative->nums, cumulative->bounds, binL2, totalDataPoints * 0.308538f);
         PS_BIN_INTERPOLATE (binH2F32, cumulative->nums, cumulative->bounds, binH2, totalDataPoints * 0.691462f);
+        PS_BIN_INTERPOLATE (binL4F32, cumulative->nums, cumulative->bounds, binL4, totalDataPoints * 0.022481f);
+        PS_BIN_INTERPOLATE (binH4F32, cumulative->nums, cumulative->bounds, binH4, totalDataPoints * 0.977519f);
 
         // report +/- 1 sigma points
@@ -874 +879 @@
                 "The exact 15.8655 and 84.1345 percent data point positions are: (%f, %f)\n",
                 binLoF32, binHiF32);
-
-        // ADD step 5: Determine SIGMA as (1/2 of) the distance between these positions.
-        // sigma = (binHiF32 - binLoF32) / 2.0;
-        sigma = (binH2F32 - binL2F32);
+        psTrace(TRACE, 5,
+                "The exact 30.8538 and 69.1462 percent data point positions are: (%f, %f)\n",
+                binL2F32, binH2F32);
+        psTrace(TRACE, 5,
+                "The exact 02.22481 and 97.7519 percent data point positions are: (%f, %f)\n",
+                binL4F32, binH4F32);
+
+        // ADD step 5: Determine SIGMA as the distance between binL2 and binH2 (+/- 0.5 sigma)
+        float sigma1 = (binH2F32 - binL2F32);
+        float sigma2 = (binHiF32 - binLoF32) / 2.0;
+        float sigma4 = (binH4F32 - binL4F32) / 4.0;
+
+        // take the smallest of the three: if we have a clump with wide outliers, sigma2 and
+        // sigma4 will be biased high; if we have a bi-modal distribution, sigma1 and sigma2
+        // will be biased high.
+        sigma = PS_MIN (sigma1, PS_MIN (sigma2, sigma4));
+
+        psTrace(TRACE, 6, "The 2x sigma is %f.\n", sigma1);
+        psTrace(TRACE, 6, "The 2x sigma is %f.\n", sigma2);
+        psTrace(TRACE, 6, "The 4x sigma is %f.\n", sigma4);
 
         psTrace(TRACE, 6, "The current sigma is %f.\n", sigma);
@@ -1698 +1719 @@
         // XXX can we calculate the binMin, binMax **before** building this histogram?
         // if the range is too absurd, adjust numBins & binSize
-        long numBins = PS_MAX (5, 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);
@@ -1846 +1867 @@
 
         // if we converge on a solution outside the range binMin - binMax, use a more conservative range
-        if (done && (guessMean > PS_BIN_MIDPOINT(histogram, binMax - 1))) {
+        float minValue = PS_BIN_MIDPOINT(histogram, binMin);
+        float maxValue = PS_BIN_MIDPOINT(histogram, binMax - 1);
+
+        if (done && ((guessMean < minValue) || (guessMean > maxValue))) {
             psTrace(TRACE, 6, "Inconsistent result, re-trying the fit\n");
 
@@ -1904 +1928 @@
             psTrace(TRACE, 6, "The symmetric stdev is %f.\n", guessStdev);
 #endif
+
+            // if we converge on a solution outside the range binMin - binMax, use a more conservative range
+            float minValueSym = PS_BIN_MIDPOINT(histogram, binS);
+            float maxValueSym = PS_BIN_MIDPOINT(histogram, binE - 1);
+
+            // saturate on min or max value
+            if (guessMean < minValueSym) {
+                guessMean = minValueSym;
+                psTrace(TRACE, 6, "The symmetric mean is out of bounds, saturating to %f.\n", guessMean);
+            }
+                
+            // saturate on min or max value
+            if (guessMean > maxValueSym) {
+                guessMean = maxValueSym;
+                psTrace(TRACE, 6, "The symmetric mean is out of bounds, saturating to %f.\n", guessMean);
+            }
 
             psFree (poly);