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

Timestamp:

Nov 22, 2006, 7:33:28 PM (20 years ago)

Author:

magnier

Message:

several important changes:

• simplified code for sampleMean and MinMax : a single copy with if statements

is only marginally slower than separate loops

• added local maskVal to two functions so external user can safely provide a 0 mask
• moved fittedStats to independent function: called multiple times if sigma changes alot
• gaussian now fits for var not sigma (avoid +/- flips)
• restrict to 2 sigma by default, but user setable

• some significant re-work of the implementation of robustStats
  • no more smoothing of the vector
  • a bit clearer choice for binsize

File:

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

trunk/psLib/src/math/psStats.c

@@ -16 +16 @@
  * use ->min and ->max (PS_STAT_USE_RANGE)
  *
- *  @version $Revision: 1.189 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-11-15 00:36:14 $
+ *  @version $Revision: 1.190 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-11-23 05:33:28 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -55 +55 @@
 /*****************************************************************************/
 #define PS_GAUSS_WIDTH 3                // The width of the Gaussian smoothing.
-// This corresponds to N in the ADD.
-#define PS_CLIPPED_NUM_ITER_LB 1
+#define PS_CLIPPED_NUM_ITER_LB 1 // This corresponds to N in the ADD.
 #define PS_CLIPPED_NUM_ITER_UB 10
 #define PS_CLIPPED_SIGMA_LB 1.0
 #define PS_CLIPPED_SIGMA_UB 10.0
 #define PS_POLY_MEDIAN_MAX_ITERATIONS 30
+#define MASK_MARK 0x80   // bit to use internally to mark data as bad
 
 #define PS_BIN_MIDPOINT(HISTOGRAM, BIN_NUM) \
@@ -104 +104 @@
 this routine sets stats->sampleMean to NAN.
  *****************************************************************************/
-static bool vectorSampleMean(const psVector* myVector,
-                             const psVector* errors,
-                             const psVector* maskVector,
-                             psMaskType maskVal,
-                             psStats* stats)
+# if (0)
+    static bool vectorSampleMean(const psVector* myVector,
+                                 const psVector* errors,
+                                 const psVector* maskVector,
+                                 psMaskType maskVal,
+                                 psStats* stats)
 {
     psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
@@ -228 +229 @@
     return true;
 }
+#endif
+
+/******************************************************************************
+vectorSampleMean(myVector, maskVector, maskVal, stats): calculates the
+mean of the input vector.  If there was a problem with the mean calculation,
+this routine sets stats->sampleMean to NAN.
+ 
+XXX : using the method below, with a single loop for various options
+      costs only a small amount and is much easier to debug. running 10000 tests
+      of 1000 point vectors for the two methods gives:
+                     separate    single 
+(mask: 0, range: 0): 0.067 sec   0.073 sec
+(mask: 1, range: 0): 0.098 sec  0.102 sec
+(mask: 0, range: 0): 0.136 sec  0.162 sec
+(mask: 1, range: 0): 0.177 sec  0.181 sec
+ *****************************************************************************/
+static bool vectorSampleMean(const psVector* myVector,
+                             const psVector* errors,
+                             const psVector* maskVector,
+                             psMaskType maskVal,
+                             psStats* stats)
+{
+    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
+
+    psF32 mean = 0.0;                   // The mean
+    long count = 0;                     // Number of points contributing to this mean
+
+    psF32 *data = myVector->data.F32;   // Dereference
+    int numData = myVector->n;          // Number of data points
+
+    psU8 *maskData = (maskVector == NULL) ? NULL : maskVector->data.U8;
+    bool useRange = stats->options & PS_STAT_USE_RANGE;
+
+    // If errors exist, we enter a slightly different loop.
+    // XXX collapse to a single loop?
+    if (!errors) {
+        for (long i = 0; i < numData; i++) {
+            // Check if the data is with the specified range
+            if (useRange && (data[i] < stats->min))
+                continue;
+            if (useRange && (data[i] > stats->max))
+                continue;
+            if (maskData && (maskData[i] & maskVal))
+                continue;
+            mean += data[i];
+            count++;
+        }
+        mean = (count > 0) ? mean / (psF32) count : NAN;
+    } else {
+        psF32 sumWeights = 0.0;         // The sum of the weights
+        psF32 *errorsData = errors->data.F32;
+        for (long i = 0; i < numData; i++) {
+            // Check if the data is with the specified range
+            if (useRange && (data[i] < stats->min))
+                continue;
+            if (useRange && (data[i] > stats->max))
+                continue;
+            if (maskData && (maskData[i] & maskVal))
+                continue;
+
+            float weight = 1.0 / PS_SQR(errorsData[i]);
+            mean += data[i] * weight;
+            sumWeights += weight;
+            count++;
+        }
+        mean = (count > 0) ? mean / sumWeights : NAN;
+    }
+
+    stats->sampleMean = mean;
+    if (isnan(mean)) {
+        psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__);
+        return false;
+    }
+
+    psTrace("psLib.math", 4, "---- %s(true) end ----\n", __func__);
+    return true;
+}
 
 /******************************************************************************
@@ -233 +311 @@
 If there was a problem with the calculation, this routine sets stats->max and stats->min to NAN.  Return the
 number of valid values in the vector (those not masked or outside the specified range.
+XXX : using the method below, with a single loop for various options
+      costs only a small amount and is much easier to debug. running 10000 tests
+      of 1000 point vectors for the two methods gives:
+                     separate    single 
+(mask: 0, range: 0):  0.101 sec  0.149 sec
+(mask: 1, range: 0):  0.125 sec  0.160 sec
+(mask: 0, range: 0):  0.179 sec  0.208 sec
+(mask: 1, range: 0):  0.200 sec  0.244 sec
  *****************************************************************************/
 static long vectorMinMax(const psVector* myVector,
@@ -244 +330 @@
     psF32 min = PS_MAX_F32;             // The calculated minimum
     psF32 *vector = myVector->data.F32; // Dereference the vector
-    psU8 *mask = NULL;                  // Dereference the mask
-    if (maskVector) {
-        mask = maskVector->data.U8;
-    }
+
     int num = myVector->n;              // Number of values
     int numValid = 0;                   // Number of valid values
 
-    // If PS_STAT_USE_RANGE is requested, then we enter a different loop.
-    if (stats->options & PS_STAT_USE_RANGE) {
-        if (mask) {
-            // Need to check range and mask
-            for (long i = 0; i < num; i++) {
-                if (!(maskVal & mask[i]) && vector[i] >= stats->min && vector[i] <= stats->max) {
-                    numValid++;
-                    if (vector[i] > max) {
-                        max = vector[i];
-                    }
-                    if (vector[i] < min) {
-                        min = vector[i];
-                    }
-                }
-            }
-        } else {
-            // Need to check range
-            for (long i = 0; i < num; i++) {
-                if (vector[i] >= stats->min && vector[i] <= stats->max) {
-                    numValid++;
-                    if (vector[i] > max) {
-                        max = vector[i];
-                    }
-                    if (vector[i] < min) {
-                        min = vector[i];
-                    }
-                }
-            }
-        }
-    } else {
-        if (mask) {
-            // Need to check mask
-            for (long i = 0; i < num; i++) {
-                if (!(maskVal & mask[i])) {
-                    numValid++;
-                    if (vector[i] > max) {
-                        max = vector[i];
-                    }
-                    if (vector[i] < min) {
-                        min = vector[i];
-                    }
-                }
-            }
-        } else {
-            // No checks
-            for (long i = 0; i < num; i++) {
-                if (vector[i] > max) {
-                    max = vector[i];
-                }
-                if (vector[i] < min) {
-                    min = vector[i];
-                }
-            }
-            numValid = num;
-        }
-    }
-
+    psU8 *maskData = (maskVector == NULL) ? NULL : maskVector->data.U8;
+    bool useRange = stats->options & PS_STAT_USE_RANGE;
+
+    for (long i = 0; i < num; i++) {
+        // Check if the data is with the specified range
+        if (useRange && (vector[i] < stats->min))
+            continue;
+        if (useRange && (vector[i] > stats->max))
+            continue;
+        if (maskData && (maskData[i] & maskVal))
+            continue;
+        numValid++;
+        max = PS_MAX (vector[i], max);
+        min = PS_MIN (vector[i], min);
+    }
 
     if (numValid == 0) {
@@ -685 +726 @@
 
 /******************************************************************************
-vectorClippedStats(myVector, errors, maskVector, maskVal, stats): calculates the
+vectorClippedStats(myVector, errors, maskInput, maskValInput, stats): calculates the
 clipped stats (mean or stdev) of the input vector.
  
 Inputs
     myVector
-    maskVector
-    maskVal
+    maskInput
+    maskValInput
     stats
 Returns
@@ -698 +739 @@
 static bool vectorClippedStats(const psVector* myVector,
                                const psVector* errors,
-                               psVector* maskVector,
-                               psMaskType maskVal,
+                               psVector* maskInput,
+                               psMaskType maskValInput,
                                psStats* stats
                               )
@@ -721 +762 @@
 
     // We copy the mask vector, to preserve the original
+    psMaskType maskVal = MASK_MARK | maskValInput;
     psVector *tmpMask = psVectorAlloc(myVector->n, PS_TYPE_U8);
     psVectorInit(tmpMask, 0);
-    if (maskVector) {
+    if (maskInput) {
         for (long i = 0; i < myVector->n; i++) {
-            if (maskVector->data.U8[i] & maskVal) {
-                tmpMask->data.U8[i] = 0xff;
+            if (maskInput->data.U8[i] & maskValInput) {
+                tmpMask->data.U8[i] = maskVal;
             }
         }
@@ -1012 +1054 @@
     psF32 x = coords->data.F32[0];
     psF32 mean = params->data.F32[0];
-    psF32 stdev = params->data.F32[1];
-
-    psF32 gauss = psGaussian(x, mean, stdev, false);
+    psF32 var = params->data.F32[1];
+    psF32 dx = (x - mean);
+
+    psF32 gauss = exp (-0.5*PS_SQR(dx)/var);
     if (deriv) {
         PS_ASSERT_VECTOR_SIZE(deriv, (long)2, NAN);
         PS_ASSERT_VECTOR_TYPE(deriv, PS_TYPE_F32, NAN);
-        psF32 tmp = (x - mean) * gauss;
-        deriv->data.F32[0] = tmp / (stdev * stdev);
-        deriv->data.F32[1] = (x - mean) * tmp / (stdev * stdev * stdev);
+        psF32 tmp = dx * gauss;
+        deriv->data.F32[0] = tmp / var;
+        deriv->data.F32[1] = tmp * dx / (var * var);
     }
 
@@ -1028 +1071 @@
 }
 
+/*
+ * vectorFittedStats requires guess for fittedMean and fittedStdev
+ * robustN50 should also be set
+ */
+static bool vectorFittedStats (const psVector* myVector,
+                               const psVector* errors,
+                               psVector* mask,
+                               psMaskType maskVal,
+                               psStats* stats)
+{
+
+    psTrace("psLib.math", 6, "The guess mean  is %f.\n", stats->fittedMean);
+    psTrace("psLib.math", 6, "The guess stdev is %f.\n", stats->fittedStdev);
+
+    psStats *statsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); // Statistics for min and max
+
+    psScalar tmpScalar;                 // Temporary scalar variable, for p_psVectorBinDisect
+    tmpScalar.type.type = PS_TYPE_F32;
+
+    // construct a histogram with (sigma/10 < binsize < sigma)
+    // set roughly so that the lowest bins have about 2 cnts
+    // Nsmallest ~ N50 / (4*dN))
+    psF32 dN = PS_MAX (1, PS_MIN (10, stats->robustN50 / 8));
+    psF32 newBinSize = stats->fittedStdev / dN;
+
+    // Determine the min/max of the vector (which prior outliers masked out)
+    int numValid = vectorMinMax(myVector, mask, maskVal, statsMinMax); // Number of values
+    float min = statsMinMax->min;
+    float max = statsMinMax->max;
+    if (numValid == 0 || isnan(min) || isnan(max)) {
+        psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
+        psFree(statsMinMax);
+        psFree(mask);
+        psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
+        return false;
+    }
+
+    // Calculate the number of bins.
+    long numBins = (max - min) / newBinSize;
+    psTrace("psLib.math", 6, "The new min/max values are (%f, %f).\n", min, max);
+    psTrace("psLib.math", 6, "The new bin size is %f.\n", newBinSize);
+    psTrace("psLib.math", 6, "The numBins is %ld\n", numBins);
+
+    psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers)
+    histogram = psVectorHistogram(histogram, myVector, errors, mask, maskVal);
+    if (psTraceGetLevel("psLib.math") >= 8) {
+        PS_VECTOR_PRINT_F32(histogram->nums);
+    }
+
+    long binMin = 0;                // Low bin to check
+    long binMax = 0;                // High bin to check
+
+    // Fit a Gaussian to the bins in the requested range about the guess mean
+    // min,max overrides clipSigma
+    psF32 maxFitSigma = 2.0;
+    if (isfinite(stats->clipSigma)) {
+        maxFitSigma = fabs(stats->clipSigma);
+    }
+    if (isfinite(stats->max)) {
+        maxFitSigma = fabs(stats->max);
+    }
+
+    psF32 minFitSigma = 2.0;
+    if (isfinite(stats->clipSigma)) {
+        minFitSigma = fabs(stats->clipSigma);
+    }
+    if (isfinite(stats->min)) {
+        minFitSigma = fabs(stats->min);
+    }
+
+    tmpScalar.data.F32 = stats->fittedMean - minFitSigma*stats->fittedStdev;
+    if (tmpScalar.data.F32 < min) {
+        binMin = 0;
+    } else {
+        binMin = p_psVectorBinDisect(histogram->bounds, &tmpScalar);
+    }
+    tmpScalar.data.F32 = stats->fittedMean + maxFitSigma*stats->fittedStdev;
+    if (tmpScalar.data.F32 > max) {
+        binMax = histogram->bounds->n - 1;
+    } else {
+        binMax = p_psVectorBinDisect(histogram->bounds, &tmpScalar);
+    }
+    if (binMin < 0 || binMax < 0) {
+        psError(PS_ERR_UNKNOWN, false, "Failed to calculate the -%f sigma : +%f sigma bins\n", minFitSigma, maxFitSigma);
+        psFree(statsMinMax);
+        psFree(histogram);
+        psFree(mask);
+        psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
+        return false;
+    }
+
+    // Generate the variables that will be used in the Gaussian fitting
+    // XXX EAM : we should test / guarantee that there are at least 3 (right?) bins
+    psVector *y = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32); // Vector of coordinates
+    psArray *x = psArrayAlloc((1 + (binMax - binMin))); // Array of ordinates
+    for (long i = binMin, j = 0; i <= binMax ; i++, j++) {
+        y->data.F32[j] = histogram->nums->data.F32[i];
+        psVector *ordinate = psVectorAlloc(1, PS_TYPE_F32); // The ordinate value
+        ordinate->data.F32[0] = PS_BIN_MIDPOINT(histogram, i);
+        x->data[j] = ordinate;
+    }
+    if (psTraceGetLevel("psLib.math") >= 8) {
+        // XXX: Print the x array somehow.
+        PS_VECTOR_PRINT_F32(y);
+    }
+    psTrace("psLib.math", 6, "The clipped numBins is %ld\n", y->n);
+
+    // Normalize y to [0.0:1.0] (since the psMinimizeLMChi2Gauss1D() functions is [0.0:1.0])
+    // XXX EAM : why not just fit the normalization as well??
+    p_psNormalizeVectorRange(y, 0.0, 1.0);
+
+    // Fit a Gaussian to the data.
+    psMinimization *minimizer = psMinimizationAlloc(100, 0.01); // The minimizer information
+    psVector *params = psVectorAlloc(2, PS_TYPE_F32); // Parameters for the Gaussian
+    // Initial guess for the mean (index 0) and var (index 1).
+    params->data.F32[0] = stats->fittedMean;
+    params->data.F32[1] = PS_SQR(stats->fittedStdev);
+    if (psTraceGetLevel("psLib.math") >= 8) {
+        PS_VECTOR_PRINT_F32(params);
+        PS_VECTOR_PRINT_F32(y);
+    }
+    if (!psMinimizeLMChi2(minimizer, NULL, params, NULL, x, y, NULL, minimizeLMChi2Gauss1D)) {
+        psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n");
+        psFree(statsMinMax);
+        psFree(x);
+        psFree(y);
+        psFree(minimizer);
+        psFree(params);
+        psFree(mask);
+        psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
+        return false;
+    }
+    if (psTraceGetLevel("psLib.math") >= 8) {
+        PS_VECTOR_PRINT_F32(params);
+    }
+
+    // The fitted mean is the Gaussian mean.
+    stats->fittedMean = params->data.F32[0];
+    psTrace("psLib.math", 6, "The fitted mean is %f.\n", stats->fittedMean);
+
+    // The fitted standard deviation
+    stats->fittedStdev = sqrt(params->data.F32[1]);
+    psTrace("psLib.math", 6, "The fitted stdev is %f.\n", stats->fittedStdev);
+
+    // Clean up after fitting
+    psFree (histogram);
+    psFree (x);
+    psFree (y);
+    psFree (minimizer);
+    psFree (params);
+    psFree (statsMinMax);
+    return true;
+}
+
+
 /******************************************************************************
-vectorRobustStats(myVector, maskVector, maskVal, stats): This is the new
+vectorRobustStats(myVector, maskInput, maskValInput, stats): This is the new
 version of the robust stats routine.
  
@@ -1049 +1247 @@
 static bool vectorRobustStats(const psVector* myVector,
                               const psVector* errors,
-                              psVector* maskVector,
-                              psMaskType maskVal,
+                              psVector* maskInput,
+                              psMaskType maskValInput,
                               psStats* stats)
 {
@@ -1061 +1259 @@
     tmpScalar.type.type = PS_TYPE_F32;
 
+    // we need to generate an internal mask and maskVal which ensures a bit is set
+    // and tested even if there is no supplied mask (and/or the maskVal is 0)
+    // XXX this would be better if we had globally defined mask values
+    psMaskType maskVal = MASK_MARK | maskValInput;
     psVector *mask = psVectorAlloc(myVector->n, PS_TYPE_MASK); // The actual mask we will use
     psVectorInit(mask, 0);
-    if (maskVector) {
+    if (maskInput) {
         for (long i = 0; i < myVector->n; i++) {
-            if (maskVector->data.U8[i] & maskVal) {
-                mask->data.U8[i] = 0xff;
+            if (maskInput->data.U8[i] & maskValInput) {
+                mask->data.U8[i] = maskVal;
             }
         }
@@ -1098 +1300 @@
         psTrace("psLib.math", 6, "Data min/max is (%.2f, %.2f)\n", min, max);
 
-        // If all data points have the same value, then we set the appropiate members of stats and return.
+        // If all data points have the same value, then we set the appropriate members of stats and return.
         if (fabs(max - min) <= FLT_EPSILON) {
             psTrace("psLib.math", 5, "All data points have the same value: %f.\n", min);
@@ -1226 +1428 @@
             return false;
         }
-
+        // XXXX we can probably just use the bin values to get sigma without interpolating...
         // ADD step 4b: Interpolate Sigma to find these two positions exactly: these are the 1sigma positions.
         #if 0
@@ -1304 +1506 @@
             for (long i = 0 ; i < myVector->n ; i++) {
                 if ((myVector->data.F32[i] < medianLo) || (myVector->data.F32[i] > medianHi)) {
-                    mask->data.U8[i] = 0xff;
+                    mask->data.U8[i] |= MASK_MARK;
                     psTrace("psLib.math", 6, "Masking element %ld is %f\n", i, myVector->data.F32[i]);
                 }
@@ -1378 +1580 @@
     psTrace("psLib.math", 6, "The robustN50 is %ld.\n", N50);
 
-
     // Fitted statistics
     if (stats->options & PS_STAT_FITTED_MEAN || stats->options & PS_STAT_FITTED_STDEV) {
-        // New algorithm for calculated mean and stdev.
-
-        // XXX: Where is this documented?
-        psF32 dN = 0.17 * N50;
-        if (dN < 1.0) {
-            dN = 1.0;
-        } else if (dN > 4.0) {
-            dN = 4.0;
-        }
-        psF32 newBinSize = sigma / dN;
-
-        // Determine the min/max of the vector (which prior outliers masked out)
-        int numValid = vectorMinMax(myVector, mask, maskVal, statsMinMax); // Number of values
-        min = statsMinMax->min;
-        max = statsMinMax->max;
-        if (numValid == 0 || isnan(min) || isnan(max)) {
-            psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
-            psFree(statsMinMax);
-            psFree(mask);
-            psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
-            return false;
-        }
-
-        // Calculate the number of bins.
-        long numBins = (max - min) / newBinSize;
-        psTrace("psLib.math", 6, "The new min/max values are (%f, %f).\n", min, max);
-        psTrace("psLib.math", 6, "The new bin size is %f.\n", newBinSize);
-        psTrace("psLib.math", 6, "The numBins is %ld\n", numBins);
-
-        psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers)
-        histogram = psVectorHistogram(histogram, myVector, errors, mask, maskVal);
-        if (psTraceGetLevel("psLib.math") >= 8) {
-            PS_VECTOR_PRINT_F32(histogram->nums);
-        }
-
-        // Smooth the resulting histogram with a Gaussian with sigma_s = 1 bin.
-        psVector *smoothed = vectorSmoothHistGaussian(histogram, newBinSize); // Smoothed histogram
-        if (psTraceGetLevel("psLib.math") >= 8) {
-            PS_VECTOR_PRINT_F32(smoothed);
-        }
-
-        // Find the bin with the peak value in the range 2 sigma of the robust histogram median.
-        long binMin = 0;                // Low bin to check
-        long binMax = 0;                // High bin to check
-        tmpScalar.data.F32 = stats->robustMedian - (2.0 * sigma);
-        if (tmpScalar.data.F32 <= histogram->bounds->data.F32[0]) {
-            binMin = 0;
-        } else {
-            binMin = p_psVectorBinDisect(histogram->bounds, &tmpScalar);
-        }
-
-        tmpScalar.data.F32 = stats->robustMedian + (2.0 + sigma);
-        if (tmpScalar.data.F32 >= histogram->bounds->data.F32[histogram->bounds->n-1]) {
-            binMax = histogram->bounds->n-1;
-        } else {
-            binMax = p_psVectorBinDisect(histogram->bounds, &tmpScalar);
-        }
-        if ((binMin < 0) || (binMax < 0)) {
-            psError(PS_ERR_UNKNOWN, false, "Failed to calculate the +/- 2.0 * sigma bins\n");
-            psFree(statsMinMax);
-            psFree(histogram);
-            psFree(mask);
-            psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
-            return false;
-        }
-
-        long binNum = binMin;           // Bin number containing the peak
-        psF32 binMaxNums = smoothed->data.F32[binNum]; // The peak value
-        for (psS32 i = binMin+1 ; i <= binMax ; i++) {
-            if (smoothed->data.F32[i] > binMaxNums) {
-                binNum = i;
-                binMaxNums = smoothed->data.F32[i];
-            }
-        }
-        psTrace("psLib.math", 6, "The peak bin is %ld, with %f data.n", binNum, binMaxNums);
-
-        // Fit a Gaussian to the bins in the range 20 sigma of the robust histogram median.
-        tmpScalar.data.F32 = stats->robustMedian - (20.0 * sigma);
-        if (tmpScalar.data.F32 < min) {
-            binMin = 0;
-        } else {
-            binMin = p_psVectorBinDisect(histogram->bounds, &tmpScalar);
-        }
-        tmpScalar.data.F32 = stats->robustMedian + (20.0 * sigma);
-        if (tmpScalar.data.F32 > max) {
-            binMax = smoothed->n - 1;
-        } else {
-            binMax = p_psVectorBinDisect(histogram->bounds, &tmpScalar);
-        }
-        if (binMin < 0 || binMax < 0) {
-            psError(PS_ERR_UNKNOWN, false, "Failed to calculate the +/- 20.0 * sigma bins\n");
-            psFree(statsMinMax);
-            psFree(histogram);
-            psFree(mask);
-            psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
-            return false;
-        }
-
-        // Generate the variables that will be used in the Gaussian fitting
-        psVector *y = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32); // Vector of coordinates
-        psArray *x = psArrayAlloc((1 + (binMax - binMin))); // Array of ordinates
-        for (long i = binMin, j = 0; i <= binMax ; i++, j++) {
-            y->data.F32[j] = smoothed->data.F32[i];
-            psVector *ordinate = psVectorAlloc(1, PS_TYPE_F32); // The ordinate value
-            ordinate->data.F32[0] = PS_BIN_MIDPOINT(histogram, i);
-            x->data[j] = ordinate;
-        }
-        if (psTraceGetLevel("psLib.math") >= 8) {
-            // XXX: Print the x array somehow.
-            PS_VECTOR_PRINT_F32(y);
-        }
-        psFree(smoothed);
-        psFree(histogram);
-
-        numValid = vectorMinMax(y, NULL, 0, statsMinMax); // Number of valid values
-        min = statsMinMax->min;
-        max = statsMinMax->max;
-        if (numValid == 0 || isnan(min) || isnan(max)) {
-            psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
-            psFree(mask);
-            psFree(histogram);
-            psFree(statsMinMax);
-            psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
-            return false;
-        }
-
-        // Normalize y to [0.0:1.0] (since the psMinimizeLMChi2Gauss1D() functions is [0.0:1.0])
-        p_psNormalizeVectorRange(y, 0.0, 1.0);
-
-        // Fit a Gaussian to the data.
-        psMinimization *minimizer = psMinimizationAlloc(100, 0.01); // The minimizer information
-        psVector *params = psVectorAlloc(2, PS_TYPE_F32); // Parameters for the Gaussian
-        // Initial guess for the mean (index 0) and standard dev (index 1).
-        params->data.F32[0] = stats->robustMedian;
-        params->data.F32[1] = sigma;
-        if (psTraceGetLevel("psLib.math") >= 8) {
-            PS_VECTOR_PRINT_F32(params);
-            PS_VECTOR_PRINT_F32(y);
-        }
-        if (!psMinimizeLMChi2(minimizer, NULL, params, NULL, x, y, NULL, minimizeLMChi2Gauss1D)) {
-            psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n");
-            psFree(statsMinMax);
-            psFree(x);
-            psFree(y);
-            psFree(minimizer);
-            psFree(params);
-            psFree(mask);
-            psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
-            return false;
-        }
-        if (psTraceGetLevel("psLib.math") >= 8) {
-            PS_VECTOR_PRINT_F32(params);
-        }
-
-        // The fitted mean is the Gaussian mean.
-        stats->fittedMean = params->data.F32[0];
-        psTrace("psLib.math", 6, "The fitted mean is %f.\n", params->data.F32[0]);
-
-        // The fitted standard deviation, SIGMA_r is determined by subtracting the smoothing scale in
-        // quadrature: SIGMA_r^2 = SIGMA^2 - sigma_s^2
-        stats->fittedStdev = sqrt(PS_SQR(params->data.F32[1]) - PS_SQR(newBinSize));
-        psTrace("psLib.math", 6, "The fitted stdev is %f.\n", stats->fittedStdev);
-
-        // Clean up after fitting
-        psFree(x);
-        psFree(y);
-        psFree(minimizer);
-        psFree(params);
+        stats->fittedStdev = stats->robustStdev;  // pass the guess sigma
+        stats->fittedMean = stats->robustMedian;  // pass the guess mean
+        vectorFittedStats (myVector, errors, mask, maskVal, stats);
+        // if there is a large swing in sigma, try a second time
+        if ((stats->fittedStdev/stats->robustStdev) < 0.75) {
+            vectorFittedStats (myVector, errors, mask, maskVal, stats);
+        }
     }
 
@@ -1914 +1954 @@
 
 /******************************************************************************
-psVectorStats(myVector, maskVector, maskVal, stats): this is the public API
+psVectorStats(in, mask, maskVal, stats): this is the public API
 function which calls the above private stats functions based on what bits
 were set in stats->options.
  
 Inputs
-    myVector
-    maskVector
+    in
+    mask
     maskVal
     stats