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

Timestamp:

Nov 1, 2004, 1:57:08 PM (22 years ago)

Author:

gusciora

Message:

Fixed the robust stats routines for stdev/mean. Fixed minimization routines
for fitting gaussians. CVS:

---

Enter Log. Lines beginning with `CVS:' are removed automatically CVS: CVS:
Committing in . CVS: CVS: Modified Files: CVS:
src/dataManip/psMinimize.c src/dataManip/psMinimize.h CVS:
src/dataManip/psStats.c src/dataManip/makedir/psMinimize.d CVS:
src/dataManip/makedir/psStats.d CVS: test/dataManip/tst_psMinimize07.c

---

File:

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

trunk/psLib/src/math/psStats.c

@@ -9 +9 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.80 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-28 22:46:57 $
+ *  @version $Revision: 1.81 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-11-01 23:57:08 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -790 +790 @@
     psFree(tmpMask);
 }
+
+void p_psNormalizeVectorF32_0(psVector* myData)
+{
+    float min = (float)HUGE;
+    float max = (float)-HUGE;
+    psS32 i = 0;
+
+    for (i = 0; i < myData->n; i++) {
+        if (myData->data.F32[i] < min) {
+            min = myData->data.F32[i];
+        }
+        if (myData->data.F32[i] > max) {
+            max = myData->data.F32[i];
+        }
+    }
+
+    //  float range = max - min;
+    //    for (i = 0; i < myData->n; i++) {
+    //        myData->data.F32[i] = (myData->data.F32[i] - min) / range;
+    //    }
+    for (i = 0; i < myData->n; i++) {
+        myData->data.F32[i] = (myData->data.F32[i] - min) / (max - min);
+    }
+}
+
+
 
 /*****************************************************************************
@@ -1162 +1188 @@
     myStdev = sqrt((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));
 
+
+    // Using the above (myMean, myStdev) as initial estimates, we fit a
+    // Gaussian to the robustHistogramVector.
+    psMinimization *min = psMinimizationAlloc(100, 0.1);
+    psVector *myParams = psVectorAlloc(2, PS_TYPE_F32);
+    psArray *myCoords = psArrayAlloc(robustHistogramVector->n);
+    psVector *y = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
+
+    p_psNormalizeVectorF32_0(robustHistogramVector);
+    for (i=0;i<robustHistogramVector->n;i++) {
+        myCoords->data[i] = (psPtr *) psVectorAlloc(2, PS_TYPE_F32);
+        ((psVector *) (myCoords->data[i]))->data.F32[0] = (float) i;
+        y->data.F32[i] = robustHistogramVector->data.F32[i];
+    }
+
+    myParams->data.F32[0] = myMean;
+    myParams->data.F32[1] = myStdev;
+    psMinimizeLMChi2(min,
+                     NULL,
+                     myParams,
+                     NULL,
+                     myCoords,
+                     y,
+                     NULL,
+                     (psMinimizeLMChi2Func) psMinimizeLMChi2Gauss1D);
+    psFree(min);
+    psFree(myParams);
+    psFree(myCoords);
+    psFree(y);
     /**************************************************************************
     Set the appropriate members in the output stats struct.
     **************************************************************************/
     if (stats->options & PS_STAT_ROBUST_MEAN) {
-        stats->robustMean = myMean;
+        if (fabs((myParams->data.F32[0] - myMean)/myMean) > 0.1) {
+            printf("WARNING: the fitted Gaussian has more than 10% error for the mean.\n");
+            printf("Using the calculated mean (calc, fit) is (%f, %f)\n", myMean, myParams->data.F32[0]);
+            stats->robustMean = myMean;
+        } else {
+            stats->robustMean = myParams->data.F32[0];
+        }
     }
 
@@ -1174 +1235 @@
 
     if (stats->options & PS_STAT_ROBUST_STDEV) {
-        stats->robustStdev = myStdev;
+        if (fabs((myParams->data.F32[1] - myStdev)/myStdev) > 0.1) {
+            printf("WARNING: the fitted Gaussian has more than 10% error for the stdev.\n");
+            printf("Using the calculated stdev (calc, fit) is (%f, %f)\n", myStdev, myParams->data.F32[1]);
+            stats->robustStdev = myStdev;
+        } else {
+            stats->robustStdev = myParams->data.F32[1];
+        }
     }