Changeset 36222 for trunk/psLib/src/math/psStats.c

Timestamp:

Oct 16, 2013, 5:49:11 PM (13 years ago)

Author:

watersc1

Message:

Partially corrected psStats code. psStats/robustMedian handles small numbers of inputs correctly. psMinimizePolyFit attempts to scale the ordinary polynomial fit data by the median and scale to make the fitting somewhat more stable. psBinomialCoeff is necessary to implement this scaling, so it is added here as a stand alone function. There is still a bias in the fitting results that I need to sort out still, but this fixes the known stats failure modes.

File:

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

trunk/psLib/src/math/psStats.c

@@ -141 +141 @@
 
 // Debug information
-#define CZW 0
+#define CZW 1
 
 /*****************************************************************************/
@@ -232 +232 @@
         }
         count++;
+
     }
     if (errors) {
@@ -424 +425 @@
     for (long i = 0; i < myVector->n; i++) {
         // Check if the data is with the specified range
+
         if (!isfinite(data[i]))
             continue;
@@ -1043 +1045 @@
 
 
-#if (CZW)
-        printf("CZW (%d): median %f sigma %f delta: %f \t %f %f %f %f %f %f %f \t %f %f %f %f %f %f %f\n",
+#if (CZW && 0)
+        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",
                iterate,
            stats->robustMedian,stats->robustStdev,
@@ -1060 +1062 @@
            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);
-        PS_VECTOR_PRINT_F32(cumulative->bounds);
-        PS_VECTOR_PRINT_F32(cumulative->nums);
+        //      PS_VECTOR_PRINT_F32(histogram->bounds);
+        //      PS_VECTOR_PRINT_F32(histogram->nums);
+        //      PS_VECTOR_PRINT_F32(cumulative->bounds);
+        //      PS_VECTOR_PRINT_F32(cumulative->nums);
 #endif
 
@@ -1276 +1278 @@
             return true;
         }
-        //      printf("BINS: %ld %f %f %f %f\n",stats->robustN50,guessStdev,binSize,min,max);
+
         // Calculate the number of bins.
         // XXX can we calculate the binMin, binMax **before** building this histogram?
@@ -1286 +1288 @@
         psTrace(TRACE, 6, "The numBins is %ld\n", numBins);
 
-        //      printf("BINS2: %ld %f %f %f %f %ld\n",stats->robustN50,guessStdev,binSize,min,max,numBins);
         
         psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers)
@@ -1385 +1386 @@
             }
             y->n = x->n = j;
-
+            
             // fit 2nd order polynomial to ln(y) = -(x-xo)^2/2sigma^2
             // XXX this fit may fail with an error for an ill-conditioned matrix (bad data)
@@ -1397 +1398 @@
                 psErrorClear();
                 COUNT_WARNING(10, 100, "Failed to fit a gaussian to the robust histogram.\n");
+
                 psFree(poly);
                 psFree(histogram);
@@ -1478 +1480 @@
             psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
             bool status = psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x);
+#if (CZW && 0)
+            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],
+                     poly->coeff[0] + poly->coeff[1] * x->data.F32[i] +
+                     poly->coeff[2] * pow(x->data.F32[i],2));
+            }
+#endif
             psFree(x);
             psFree(y);
@@ -1493 +1503 @@
             fullfitStdev = sqrt(-0.5/poly->coeff[2]);
             fullfitMean = poly->coeff[1]*PS_SQR(fullfitStdev);
+
 #ifndef PS_NO_TRACE
             psTrace(TRACE, 6, "Parabolic Symmetric fit results: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]);
@@ -1515 +1526 @@
             }
 
+            
             psFree (poly);
         }
@@ -1537 +1549 @@
             done = true;
         }
-#if (CZW) 
-        printf("CZW IN FITTED? low %f %f full %f %f robust %f %f final %f %f\n",
-               lowfitMean,lowfitStdev,fullfitMean,fullfitStdev,stats->robustMedian,stats->robustStdev,
+
+        
+#if (CZW && 0)
+        printf("CZW IN FITTED: iter   %d %f \n"
+               "               low    %f %f \n"
+               "               full   %f %f \n"
+               "               robust %f %f \n"
+               "               final  %f %f\n",
+               iteration,minValueSym,
+               lowfitMean,lowfitStdev,
+               fullfitMean,fullfitStdev,
+               stats->robustMedian,stats->robustStdev,
                guessMean,guessStdev);
 #endif
@@ -2423 +2444 @@
         psTrace(TRACE, 6, "y data is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]);
 
-#if (CZW)
-        printf("  polyin: %f %f %f == %f %f %f\n",
-               x->data.F64[0], x->data.F64[1], x->data.F64[2],
-               y->data.F64[0], y->data.F64[1], y->data.F64[2]);
-        printf("  rawpolyin: %f %f %f == %f %f %f\n",
-               xVec->data.F32[binNum - 1], xVec->data.F32[binNum], xVec->data.F32[binNum + 1],
-               y->data.F64[0], y->data.F64[1], y->data.F64[2]);
-#endif
 
         // Ensure that the y value lies within range of the y values.
@@ -2469 +2482 @@
                 (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[1]),
                 (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[2]));
-#if (CZW)
-        printf("  poly: %f %f %f fit: %f %f %f\n",
-               myPoly->coeff[0],myPoly->coeff[1],myPoly->coeff[2],
-               (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[0]),
-               (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[1]),
-               (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[2]));
-#endif
 
         psTrace(TRACE, 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal);
@@ -2498 +2504 @@
         //        tmpFloat = xVec->data.F32[fitBin] + dY * dX;
         tmpFloat = binValue;
-                
-#if (CZW)
-        printf("   internal median: %f %f\n",tmpFloat,binValue);
-#endif
         
     } else {
@@ -2556 +2558 @@
     )
 {
+
 # if (1)
-    double Sx = (yVec->data.F32[binNum - 2] + yVec->data.F32[binNum - 1] +
-                 yVec->data.F32[binNum - 0] +
-                 yVec->data.F32[binNum + 1] + yVec->data.F32[binNum + 2]);
-
-    double YM = 0.5*(xVec->data.F32[binNum - 2] + xVec->data.F32[binNum - 1]);
-    double Ym = 0.5*(xVec->data.F32[binNum - 1] + xVec->data.F32[binNum - 0]);
-    double Yo = 0.5*(xVec->data.F32[binNum - 0] + xVec->data.F32[binNum + 1]);
-    double Yp = 0.5*(xVec->data.F32[binNum + 1] + xVec->data.F32[binNum + 2]);
-    double YP = 0.5*(xVec->data.F32[binNum + 2] + xVec->data.F32[binNum + 3]);
-
-    double Sy = YM + Ym + Yo + Yp + YP;
-
-    double Sxx = (PS_SQR(yVec->data.F32[binNum - 2]) + PS_SQR(yVec->data.F32[binNum - 1]) +
-                  PS_SQR(yVec->data.F32[binNum - 0]) +
-                  PS_SQR(yVec->data.F32[binNum + 1]) + PS_SQR(yVec->data.F32[binNum + 2]));
-
-    double Sxy = (yVec->data.F32[binNum - 2] * YM +
-                  yVec->data.F32[binNum - 1] * Ym +
-                  yVec->data.F32[binNum - 0] * Yo +
-                  yVec->data.F32[binNum + 1] * Yp +
-                  yVec->data.F32[binNum + 2] * YP);
-
-    double Det = 5*Sxx - Sx*Sx;
-    if (Det == 0.0) return NAN;
-
-    double C0 = (Sy*Sxx - Sx*Sxy) / Det;
-    double C1 = (Sxy*5 - Sx*Sy) / Det;
-
-    double value = C0 + yVal*C1;
-
-    return value;
+# define HALF_SIZE 2
+  double Sx = 0.0;
+
+  double Sy = 0.0;
+  double Sxx = 0.0;
+  double Sxy = 0.0;
+  double deltaY = 0.0;
+  int N = 0;
+
+  for (int u = binNum - HALF_SIZE; u <= binNum + HALF_SIZE; u++) {
+    if ((u >= 0)&&(u < yVec->n)) {
+      if (u+1 < xVec->n) {
+        Sx += yVec->data.F32[u];
+        Sxx += PS_SQR(yVec->data.F32[u]);
+
+        deltaY = 0.5 * (xVec->data.F32[u] + xVec->data.F32[u+1]);
+        Sy += deltaY;
+        Sxy += yVec->data.F32[u] * deltaY;
+        N += 1;
+      }
+    }
+  }
+  double Det = N * Sxx - Sx * Sx;
+  if (Det == 0.0) return NAN;
+  if (N == 0) return NAN;
+
+  double C0 = (Sy*Sxx - Sx*Sxy) / Det;
+  double C1 = (Sxy*N - Sx*Sy) / Det;
+  
+  double value = C0 + yVal*C1;
+  return value;
+  
+  
 # else
     psTrace(TRACE, 5, "binNum, yVal is (%d, %f)\n", binNum, yVal);
@@ -2682 +2687 @@
         tmpFloat = binValue;
                 
-#if (CZW)
-        printf("   internal median: %f %f\n",tmpFloat,binValue);
-#endif
         
     } else {