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

Timestamp:

Jun 23, 2004, 10:45:16 AM (22 years ago)

Author:

gusciora

Message:

Coded recent ADD changes into the robust stats. Still waiting on
minimization routines to be coded.

File:

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

trunk/psLib/src/math/psStats.c

@@ -847 +847 @@
 
 /******************************************************************************
+p_psVectorSampleStdev(myVector, maskVector, maskVal, stats): calculates the
+stdev of the input vector.
+Inputs
+    myVector
+    maskVector
+    maskVal
+    stats
+Returns
+    NULL
+ 
+NOTE: the mean is always calculated exactly.  Robust means are never
+calculated in this routine.
+ *****************************************************************************/
+void p_psVectorSampleStdev(const psVector *restrict myVector,
+                           const psVector *restrict maskVector,
+                           unsigned int maskVal,
+                           psStats *stats)
+{
+    int i = 0;                                  // Loop index variable
+    int countInt = 0;                           // # of data points being used
+    float countFloat = 0.0;                     // # of data points being used
+    float mean = 0.0;                           // The mean
+    float diff = 0.0;                           // Used in calculating stdev
+    float sumSquares = 0.0;                     // temporary variable
+    float sumDiffs = 0.0;                       // temporary variable
+    float rangeMin = 0.0;                       // Exclude data below this
+    float rangeMax = 0.0;                       // Exclude date above this
+
+    // This procedure requires the mean.  If it has not been already
+    // calculated, then call p_psVectorSampleMean()
+    if (0 != isnan(stats->sampleMean)) {
+        p_psVectorSampleMean(myVector, maskVector, maskVal, stats);
+    }
+    mean = stats->sampleMean;
+
+    if (stats->options & PS_STAT_USE_RANGE) {
+        if (maskVector != NULL) {
+            for (i=0;i<myVector->n;i++) {
+                if (!(maskVal & maskVector->data.U8[i]) &&
+                        (rangeMin <= myVector->data.F32[i]) &&
+                        (myVector->data.F32[i] <= rangeMax)) {
+                    diff = myVector->data.F32[i] - mean;
+                    sumSquares+= (diff * diff);
+                    sumDiffs+= diff;
+                    countInt++;
+                }
+            }
+        } else {
+            for (i=0;i<myVector->n;i++) {
+                if ((rangeMin <= myVector->data.F32[i]) &&
+                        (myVector->data.F32[i] <= rangeMax)) {
+                    diff = myVector->data.F32[i] - mean;
+                    sumSquares+= (diff * diff);
+                    sumDiffs+= diff;
+                    countInt++;
+                }
+            }
+            countInt = myVector->n;
+        }
+    } else {
+        if (maskVector != NULL) {
+            for (i=0;i<myVector->n;i++) {
+                if (!(maskVal & maskVector->data.U8[i])) {
+                    diff = myVector->data.F32[i] - mean;
+                    sumSquares+= (diff * diff);
+                    sumDiffs+= diff;
+                    countInt++;
+                }
+            }
+        } else {
+            for (i=0;i<myVector->n;i++) {
+                diff = myVector->data.F32[i] - mean;
+                sumSquares+= (diff * diff);
+                sumDiffs+= diff;
+                countInt++;
+            }
+            countInt = myVector->n;
+        }
+    }
+    countFloat = (float) countInt;
+
+    #ifdef DARWIN
+
+    stats->sampleStdev = (float) sqrt( (sumSquares-(sumDiffs *
+                                        sumDiffs/countFloat))/ (countFloat-1));
+    #else
+
+    stats->sampleStdev = sqrtf( (sumSquares-(sumDiffs *
+                                 sumDiffs/countFloat))/ (countFloat-1));
+    #endif
+}
+
+/******************************************************************************
+p_psVectorClippedStats(myVector, maskVector, maskVal, stats): calculates the
+clipped stats (mean or stdev) of the input vector.
+ 
+Inputs
+    myVector
+    maskVector
+    maskVal
+    stats
+Returns
+    NULL
+ *****************************************************************************/
+void p_psVectorClippedStats(const psVector *restrict myVector,
+                            const psVector *restrict maskVector,
+                            unsigned int maskVal,
+                            psStats *stats)
+{
+    int i = 0;                                  // Loop index variable
+    int j = 0;                                  // Loop index variable
+    float clippedMean = 0.0;                    // self-explanatory
+    float clippedStdev = 0.0;                   // self-explanatory
+    float oldStanMean = 0.0;                    // Temporary variable
+    float oldStanStdev = 0.0;                   // Temporary variable
+    psVector *tmpMask = NULL;                   // Temporary vector
+
+    // Endure that stats->clipIter is within the proper range.
+    if (!((CLIPPED_NUM_ITER_LB <= stats->clipIter ) &&
+            (stats->clipIter <= CLIPPED_NUM_ITER_UB))) {
+        psAbort(__func__, "Unallowed value for clipIter (%d).\n",
+                stats->clipIter);
+    }
+
+    // Endure that stats->clipSigma is within the proper range.
+    if (!((CLIPPED_SIGMA_LB <= stats->clipSigma ) &&
+            (stats->clipSigma <= CLIPPED_SIGMA_UB))) {
+        psAbort(__func__, "Unallowed value for clipSigma (%f).\n",
+                stats->clipSigma);
+    }
+
+    // We allocate a temporary mask vector since during the iterative
+    // steps that follow, we will be masking off additional data points.
+    // However, we do no want to modify the original mask vector.
+    tmpMask = psVectorAlloc(myVector->n, PS_TYPE_U8);
+    tmpMask->n = myVector->n;
+
+    // If we were called with a mask vector, then initialize the temporary
+    // mask vector with those values.
+    if (maskVector != NULL) {
+        for (i=0;i<tmpMask->n;i++) {
+            tmpMask->data.U8[i] = maskVector->data.U8[i];
+        }
+    }
+
+    // 1. Compute the sample median.
+    // NOTE: This seems odd.  Verify with IfA that we want to calculate the
+    // median here, not the mean.
+    p_psVectorSampleMedian(myVector, maskVector, maskVal, stats);
+
+    // 2. Compute the sample standard deviation.
+    p_psVectorSampleStdev(myVector, maskVector, maskVal, stats);
+
+    // 3. Use the sample median as the first estimator of the mean X.
+    clippedMean = stats->sampleMean;
+
+    // 4. Use the sample stdev as the first estimator of the mean stdev.
+    clippedStdev = stats->sampleStdev;
+
+    // Must save the old sampleMean and sampleStdev since the following code
+    // block overwrites them.
+    oldStanMean = stats->sampleMean;
+    oldStanStdev = stats->sampleStdev;
+
+    // 5. Repeat N times:
+    for (i=0;i<stats->clipIter;i++) {
+        for (j=0;j<myVector->n;j++) {
+            // a) Exclude all values x_i for which |x_i - x| > K * stdev
+            if (fabs(myVector->data.F32[j] - clippedMean) >
+                    (stats->clipSigma * clippedStdev)) {
+                tmpMask->data.U8[i] = 0xff;
+            }
+            // b) compute new mean and stdev
+            p_psVectorSampleMedian(myVector, tmpMask, maskVal, stats);
+            p_psVectorSampleStdev(myVector, tmpMask, maskVal, stats);
+
+            // c) Use the new mean for x
+            clippedMean = stats->sampleMean;
+
+            // d) Use the new stdev for stdev
+            clippedStdev = stats->sampleStdev;
+        }
+    }
+    stats->sampleMean = oldStanMean;
+    stats->sampleStdev= oldStanStdev;
+
+    // 7. The last calcuated value of x is the cliped mean.
+    if (stats->options & PS_STAT_CLIPPED_MEAN) {
+        stats->clippedMean = clippedMean;
+    }
+
+    // 8. The last calcuated value of stdev is the cliped stdev.
+    if (stats->options & PS_STAT_CLIPPED_STDEV) {
+        stats->clippedStdev = clippedStdev;
+    }
+
+    psVectorFree(tmpMask);
+}
+
+
+/******************************************************************************
 p_psVectorRobustStats(myVector, maskVector, maskVal, stats): this procedure
 calculates a variety of robust stat measures:
@@ -884 +1085 @@
     float dL = 0.0;
     int numBins = 0;
+    psStats *tmpStats = psStatsAlloc(PS_STAT_CLIPPED_STDEV|PS_STAT_CLIPPED_MEAN);
 
     // NOTE: The SDRS states that the sample quartiles must be used to
@@ -890 +1092 @@
     // size of the data set.  We should consult with IfA to ensure that this
     // is really required.
-
-    if (isnan(stats->sampleUQ) ||
+    /*
+        if (isnan(stats->sampleUQ) ||
             isnan(stats->sampleLQ)) {
-        stats->options = stats->options | PS_STAT_SAMPLE_QUARTILE;
-        p_psVectorSampleQuartiles(myVector,
-                                  maskVector,
-                                  maskVal,
-                                  stats);
-    }
-
+            stats->options = stats->options | PS_STAT_SAMPLE_QUARTILE;
+            p_psVectorSampleQuartiles(myVector,
+                                      maskVector,
+                                      maskVal,
+                                      stats);
+        }
+    */
     // Compute the initial bin size of the robust histogram.
-    sigmaE = (stats->sampleUQ - stats->sampleLQ) / 1.34f;
-    binSize = sigmaE / 10.0f;
+    p_psVectorClippedStats(myVector, maskVector, maskVal, tmpStats);
+    binSize = stats->clippedStdev / 10.0f;
+
+    // If stats->clippedStdev == 0.0, then all data elements have the same
+    // value.  Therefore, we can set the appropiate results and return.
+    if (fabs(binSize) <= FLT_EPSILON) {
+        if (stats->options & PS_STAT_ROBUST_MEAN) {
+            stats->robustMean = stats->clippedMean;
+        }
+        if  (stats->options & PS_STAT_ROBUST_MEDIAN) {
+            stats->robustMedian = stats->clippedMean;
+        }
+        if  (stats->options & PS_STAT_ROBUST_MODE) {
+            stats->robustMode = stats->clippedMean;
+        }
+        if  (stats->options & PS_STAT_ROBUST_STDEV) {
+            stats->robustStdev = 0.0;
+        }
+        if  (stats->options & PS_STAT_ROBUST_QUARTILE) {
+            stats->robustUQ = stats->clippedMean;
+            stats->robustLQ = stats->clippedMean;
+        }
+        psStatsFree(tmpStats);
+        psHistogramFree(robustHistogram);
+        return;
+    }
 
     // Detemine minimum and maximum values in the data vector.
@@ -912 +1138 @@
     }
 
-    // Create the histogram structure (yes, 2 is necessary, not 1).
+    // Create the histogram structure (yes, 2 is necessary, not 1).  Also,
+    // if we get here, we know that binSize != 0.0.
     numBins = 2 + (int) ((stats->max - stats->min) / binSize);
 
@@ -973 +1200 @@
     stats->robustNfit = 0.0;
     stats->robustN50 = 0.0;
+    psStatsFree(tmpStats);
     psHistogramFree(robustHistogram);
-}
-
-/******************************************************************************
-p_psVectorSampleStdev(myVector, maskVector, maskVal, stats): calculates the
-stdev of the input vector.
-Inputs
-    myVector
-    maskVector
-    maskVal
-    stats
-Returns
-    NULL
- 
-NOTE: the mean is always calculated exactly.  Robust means are never
-calculated in this routine.
- *****************************************************************************/
-void p_psVectorSampleStdev(const psVector *restrict myVector,
-                           const psVector *restrict maskVector,
-                           unsigned int maskVal,
-                           psStats *stats)
-{
-    int i = 0;                                  // Loop index variable
-    int countInt = 0;                           // # of data points being used
-    float countFloat = 0.0;                     // # of data points being used
-    float mean = 0.0;                           // The mean
-    float diff = 0.0;                           // Used in calculating stdev
-    float sumSquares = 0.0;                     // temporary variable
-    float sumDiffs = 0.0;                       // temporary variable
-    float rangeMin = 0.0;                       // Exclude data below this
-    float rangeMax = 0.0;                       // Exclude date above this
-
-    // This procedure requires the mean.  If it has not been already
-    // calculated, then call p_psVectorSampleMean()
-    if (0 != isnan(stats->sampleMean)) {
-        p_psVectorSampleMean(myVector, maskVector, maskVal, stats);
-    }
-    mean = stats->sampleMean;
-
-    if (stats->options & PS_STAT_USE_RANGE) {
-        if (maskVector != NULL) {
-            for (i=0;i<myVector->n;i++) {
-                if (!(maskVal & maskVector->data.U8[i]) &&
-                        (rangeMin <= myVector->data.F32[i]) &&
-                        (myVector->data.F32[i] <= rangeMax)) {
-                    diff = myVector->data.F32[i] - mean;
-                    sumSquares+= (diff * diff);
-                    sumDiffs+= diff;
-                    countInt++;
-                }
-            }
-        } else {
-            for (i=0;i<myVector->n;i++) {
-                if ((rangeMin <= myVector->data.F32[i]) &&
-                        (myVector->data.F32[i] <= rangeMax)) {
-                    diff = myVector->data.F32[i] - mean;
-                    sumSquares+= (diff * diff);
-                    sumDiffs+= diff;
-                    countInt++;
-                }
-            }
-            countInt = myVector->n;
-        }
-    } else {
-        if (maskVector != NULL) {
-            for (i=0;i<myVector->n;i++) {
-                if (!(maskVal & maskVector->data.U8[i])) {
-                    diff = myVector->data.F32[i] - mean;
-                    sumSquares+= (diff * diff);
-                    sumDiffs+= diff;
-                    countInt++;
-                }
-            }
-        } else {
-            for (i=0;i<myVector->n;i++) {
-                diff = myVector->data.F32[i] - mean;
-                sumSquares+= (diff * diff);
-                sumDiffs+= diff;
-                countInt++;
-            }
-            countInt = myVector->n;
-        }
-    }
-    countFloat = (float) countInt;
-
-    #ifdef DARWIN
-
-    stats->sampleStdev = (float) sqrt( (sumSquares-(sumDiffs *
-                                        sumDiffs/countFloat))/ (countFloat-1));
-    #else
-
-    stats->sampleStdev = sqrtf( (sumSquares-(sumDiffs *
-                                 sumDiffs/countFloat))/ (countFloat-1));
-    #endif
-}
-
-/******************************************************************************
-p_psVectorClippedStats(myVector, maskVector, maskVal, stats): calculates the
-clipped stats (mean or stdev) of the input vector.
- 
-Inputs
-    myVector
-    maskVector
-    maskVal
-    stats
-Returns
-    NULL
- *****************************************************************************/
-void p_psVectorClippedStats(const psVector *restrict myVector,
-                            const psVector *restrict maskVector,
-                            unsigned int maskVal,
-                            psStats *stats)
-{
-    int i = 0;                                  // Loop index variable
-    int j = 0;                                  // Loop index variable
-    float clippedMean = 0.0;                    // self-explanatory
-    float clippedStdev = 0.0;                   // self-explanatory
-    float oldStanMean = 0.0;                    // Temporary variable
-    float oldStanStdev = 0.0;                   // Temporary variable
-    psVector *tmpMask = NULL;                   // Temporary vector
-
-    // Endure that stats->clipIter is within the proper range.
-    if (!((CLIPPED_NUM_ITER_LB <= stats->clipIter ) &&
-            (stats->clipIter <= CLIPPED_NUM_ITER_UB))) {
-        psAbort(__func__, "Unallowed value for clipIter (%d).\n",
-                stats->clipIter);
-    }
-
-    // Endure that stats->clipSigma is within the proper range.
-    if (!((CLIPPED_SIGMA_LB <= stats->clipSigma ) &&
-            (stats->clipSigma <= CLIPPED_SIGMA_UB))) {
-        psAbort(__func__, "Unallowed value for clipSigma (%f).\n",
-                stats->clipSigma);
-    }
-
-    // We allocate a temporary mask vector since during the iterative
-    // steps that follow, we will be masking off additional data points.
-    // However, we do no want to modify the original mask vector.
-    tmpMask = psVectorAlloc(myVector->n, PS_TYPE_U8);
-    tmpMask->n = myVector->n;
-
-    // If we were called with a mask vector, then initialize the temporary
-    // mask vector with those values.
-    if (maskVector != NULL) {
-        for (i=0;i<tmpMask->n;i++) {
-            tmpMask->data.U8[i] = maskVector->data.U8[i];
-        }
-    }
-
-    // 1. Compute the sample median.
-    // NOTE: This seems odd.  Verify with IfA that we want to calculate the
-    // median here, not the mean.
-    p_psVectorSampleMedian(myVector, maskVector, maskVal, stats);
-
-    // 2. Compute the sample standard deviation.
-    p_psVectorSampleStdev(myVector, maskVector, maskVal, stats);
-
-    // 3. Use the sample median as the first estimator of the mean X.
-    clippedMean = stats->sampleMean;
-
-    // 4. Use the sample stdev as the first estimator of the mean stdev.
-    clippedStdev = stats->sampleStdev;
-
-    // Must save the old sampleMean and sampleStdev since the following code
-    // block overwrites them.
-    oldStanMean = stats->sampleMean;
-    oldStanStdev = stats->sampleStdev;
-
-    // 5. Repeat N times:
-    for (i=0;i<stats->clipIter;i++) {
-        for (j=0;j<myVector->n;j++) {
-            // a) Exclude all values x_i for which |x_i - x| > K * stdev
-            if (fabs(myVector->data.F32[j] - clippedMean) >
-                    (stats->clipSigma * clippedStdev)) {
-                tmpMask->data.U8[i] = 0xff;
-            }
-            // b) compute new mean and stdev
-            p_psVectorSampleMedian(myVector, tmpMask, maskVal, stats);
-            p_psVectorSampleStdev(myVector, tmpMask, maskVal, stats);
-
-            // c) Use the new mean for x
-            clippedMean = stats->sampleMean;
-
-            // d) Use the new stdev for stdev
-            clippedStdev = stats->sampleStdev;
-        }
-    }
-    stats->sampleMean = oldStanMean;
-    stats->sampleStdev= oldStanStdev;
-
-    // 7. The last calcuated value of x is the cliped mean.
-    if (stats->options & PS_STAT_CLIPPED_MEAN) {
-        stats->clippedMean = clippedMean;
-    }
-
-    // 8. The last calcuated value of stdev is the cliped stdev.
-    if (stats->options & PS_STAT_CLIPPED_STDEV) {
-        stats->clippedStdev = clippedStdev;
-    }
-
-    psVectorFree(tmpMask);
 }