Changeset 10550

Timestamp:

Dec 8, 2006, 1:38:54 AM (19 years ago)

Author:

magnier

Message:

• moved psHistogram functions to psHistogram.c,h (pslib_strict.h, MAkefile.am, psMemory.c)
• changed psVectorStats return value to bool (psImageGeomManip.c, psImagePixelExtract.c, psImageStats.c, psMinPoly.c,
• added FITTED_MEAN,STDEV_V2 (quadratic fit to log-peak)
• added stats->results
• set stats->results when stats are measured
• return stats calculated even if not requested
• moved vectorFittedStats call out of vectorRobustStats to psVectorStats
• internal calls to private stats function XXX now tests stats->results for the corresponding flag, NOT if value of stats->XXX == NAN.
• dropped unused function NonEmpty (commented out)
• rolled up vectorSampleMedian with if tests in the loop
• rolled up vectorSampleStdev with if tests in the loop
• call psError on failure for any internal stats call (the user should call psErrorClear if they ignore errors from psVectorStats)

Location:

trunk/psLib

Files:

• src/imageops/psImageBackground.c (modified) (2 diffs)
• src/math/psHistogram.c (added)
• src/math/psHistogram.h (added)
• src/math/psStats.c (modified) (47 diffs)
• src/math/psStats.h (modified) (7 diffs)
• test/math/tap_psStatsTiming.c (modified) (10 diffs)
• test/math/tap_psStatsTiming.txt (added)
• test/math/tap_psStats_Sample_01.c (modified) (3 diffs)

trunk/psLib/src/imageops/psImageBackground.c

@@ -73 +73 @@
         if (!psVectorSort(values, values)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to sort values.\n");
-            psFree(stats);
             psFree(values);
-            return NULL;
+            return false;
         }
 
@@ -91 +90 @@
     } else {
         // XXX leave this as a psphot user option (passed in as part of stats?)
-        if (stats->options & PS_STAT_FITTED_MEAN) {
+        if ((stats->options & PS_STAT_FITTED_MEAN) || (stats->options & PS_STAT_FITTED_MEAN_V2)) {
             stats->clipSigma = 1.0;
         }
-        if (psVectorStats (stats, values, NULL, NULL, 0) == NULL) {
+        if (!psVectorStats (stats, values, NULL, NULL, 0)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to measure statistics for image background "
                     "(%dx%d, (row0,col0) = (%d,%d)",
                     image->numRows, image->numCols, image->row0, image->col0);
-            psFree(stats);
             psFree(values);
-            return NULL;
+            return false;
         }
     }
 
     psFree(values);
-    return stats;
+    return true;
 }

trunk/psLib/src/math/psStats.c

@@ -7 +7 @@
  *  on those data structures.
  *
- *  @author GLG, MHPCC
- *
- *  XXX: The following stats members are never used, or set in this code.
- *      stats->clippedNvalues
+ *  @author GLG (MHPCC), EAM (IfA)
  *
  * XXX: Must do
@@ -16 +13 @@
  * use ->min and ->max (PS_STAT_USE_RANGE)
  *
- *  @version $Revision: 1.193 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-12-05 20:05:57 $
+ *  @version $Revision: 1.194 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-12-08 11:38:54 $
  *
- *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
+ *  Copyright 2006 IfA, University of Hawaii
  */
 
@@ -35 +32 @@
 #include "psMemory.h"
 #include "psAbort.h"
-#include "psImage.h"
+//#include "psImage.h"
 #include "psVector.h"
 #include "psTrace.h"
 #include "psLogMsg.h"
 #include "psError.h"
+#include "psHistogram.h"
 #include "psStats.h"
 #include "psMinimizeLMM.h"
@@ -64 +62 @@
 #define PS_BIN_MIDPOINT(HISTOGRAM, BIN_NUM) \
 (0.5 * (HISTOGRAM->bounds->data.F32[(BIN_NUM)] + HISTOGRAM->bounds->data.F32[(BIN_NUM)+1]))
+
 /*****************************************************************************/
 /* TYPE DEFINITIONS                                                          */
@@ -87 +86 @@
 NOTE: it is assumed that any call to these statistical functions will have
 been preceded by a call to the psVectorStats() function.  Various sanity tests
-will only be performed in psVectorStats().  Should we perform the sanity
-checks in each routine anyway?
+will only be performed in psVectorStats().
  
-XXX: For many of these private stats routines, what should be done if there
-are no acceptable elements in the input vector (if no elements lie within
-range, or there are no unmasked elements, or the input vector is NULL)?
-Currently we set the value to NAN.
- 
-XXX: Optimization: many routines have an "empty" boolean variable which keeps
-track of whether or not the vector has any valid elements.  This code can
-possibly be optimized away.
- *****************************************************************************/
-/******************************************************************************
-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.
- *****************************************************************************/
-# 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__);
-
-    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
-
-    // If PS_STAT_USE_RANGE is requested, then we enter a slightly different loop.
-    if (!errors) {
-        if (stats->options & PS_STAT_USE_RANGE) {
-            if (maskVector) {
-                // No errors, check range, check mask
-                psU8 *maskData = maskVector->data.U8;
-                for (long i = 0; i < numData; i++) {
-                    // Check if the data is with the specified range
-                    if (!(maskVal & maskData[i]) && stats->min <= data[i] && data[i] <= stats->max) {
-                        mean += data[i];
-                        count++;
-                    }
-                }
-            } else {
-                // No errors, check range, no mask
-                for (long i = 0; i < numData; i++) {
-                    if (stats->min <= data[i] && data[i] <= stats->max) {
-                        mean += data[i];
-                        count++;
-                    }
-                }
-            }
-        } else {
-            if (maskVector) {
-                // No errors, no range check, check mask
-                psU8 *maskData = maskVector->data.U8;
-                for (long i = 0; i < numData; i++) {
-                    if (!(maskVal & maskData[i])) {
-                        mean += data[i];
-                        count++;
-                    }
-                }
-
-            } else {
-                // No errors, no range check, no mask
-                for (long i = 0; i < numData; i++) {
-                    mean += data[i];
-                }
-                count = numData;
-            }
-        }
-
-        if (count > 0) {
-            mean /= (psF32)count;
-        } else {
-            mean = NAN;
-        }
-
-    } else {
-        psF32 sumWeights = 0.0;         // The sum of the weights
-        psF32 *errorsData = errors->data.F32;
-        if (stats->options & PS_STAT_USE_RANGE) {
-            if (maskVector) {
-                psU8 *maskData = maskVector->data.U8;
-                for (long i = 0; i < numData; i++) {
-                    // Check if the data is with the specified range
-                    if (!(maskVal & maskData[i]) && stats->min <= data[i] && data[i] <= stats->max) {
-                        float weight = 1.0 / PS_SQR(errorsData[i]);
-                        mean += data[i] * weight;
-                        sumWeights += weight;
-                        count++;
-                    }
-                }
-            } else {
-                for (long i = 0; i < myVector->n; i++) {
-                    if (stats->min <= data[i] && data[i] <= stats->max) {
-                        float weight = 1.0 / PS_SQR(errorsData[i]);
-                        mean += data[i] * weight;
-                        sumWeights += weight;
-                        count++;
-                    }
-                }
-            }
-        } else {
-            if (maskVector) {
-                psU8 *maskData = maskVector->data.U8;
-                for (long i = 0; i < myVector->n; i++) {
-                    if (!(maskVal & maskData[i])) {
-                        float weight = 1.0 / PS_SQR(errorsData[i]);
-                        mean += data[i] * weight;
-                        sumWeights += weight;
-                        count++;
-                    }
-                }
-            } else {
-                for (long i = 0; i < myVector->n; i++) {
-                    float weight = 1.0 / PS_SQR(errorsData[i]);
-                    mean += data[i] * weight;
-                    sumWeights += weight;
-                }
-                count = myVector->n;
-            }
-        }
-
-        if (count > 0) {
-            mean /= sumWeights;
-        } else {
-            mean = 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;
-}
-#endif
+For many of these private stats routines, it is possible that there are no acceptable elements
+in the input vector (if no elements lie within range, or there are no unmasked elements, or the
+input vector is NULL).  In such cases, we set the desired value to NAN and return an error.
+The calling function may clear the error.
+*****************************************************************************/
+
+// static prototypes:
+static psF32 minimizeLMChi2Gauss1D(psVector *deriv, const psVector *params, const psVector *coords);
+static psF32 fitQuadraticSearchForYThenReturnX(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal);
 
 /******************************************************************************
@@ -236 +103 @@
 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
- *****************************************************************************/
+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          w/errors
+(mask: 0, range: 0): 0.067 sec  0.073 sec (10%)  0.072 sec
+(mask: 1, range: 0): 0.098 sec  0.102 sec ( 4%)  0.119 sec
+(mask: 0, range: 0): 0.136 sec  0.162 sec (20%)  0.170 sec
+(mask: 1, range: 0): 0.177 sec  0.181 sec ( 2%)  0.198 sec
+ 
+(effect of errors not tested)
+ 
+To optmize this, use a macro and ifdef in or out the three states (errors, mask, range)
+*****************************************************************************/
 static bool vectorSampleMean(const psVector* myVector,
                              const psVector* errors,
@@ -253 +124 @@
     psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
 
+    long count = 0;                     // Number of points contributing to this mean
     psF32 mean = 0.0;                   // The mean
-    long count = 0;                     // Number of points contributing to this mean
+    psF32 weight;
 
     psF32 *data = myVector->data.F32;   // Dereference
@@ -262 +134 @@
     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]);
+    psF32 sumWeights = 0.0;  // The sum of the weights
+    psF32 *errorsData = (errors == NULL) ? NULL : 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;
+        if (errors) {
+            weight = (errorsData[i] == 0) ? 0.0 : PS_SQR(errorsData[i]);
             mean += data[i] * weight;
             sumWeights += weight;
-            count++;
-        }
+        } else {
+            mean += data[i];
+        }
+        count++;
+    }
+    if (errors) {
         mean = (count > 0) ? mean / sumWeights : NAN;
-    }
-
+    } else {
+        mean = (count > 0) ? mean / count : NAN;
+    }
     stats->sampleMean = mean;
     if (isnan(mean)) {
+        // XXX raise an error here?
         psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__);
         return false;
     }
 
+    stats->results |= PS_STAT_SAMPLE_MEAN;
     psTrace("psLib.math", 4, "---- %s(true) end ----\n", __func__);
     return true;
@@ -319 +183 @@
 (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,
                          const psVector* maskVector,
@@ -350 +214 @@
     }
 
+    // XXX save numValid in psStats?
     if (numValid == 0) {
         stats->max = NAN;
@@ -356 +221 @@
         stats->max = max;
         stats->min = min;
+        stats->results |= PS_STAT_MIN;
+        stats->results |= PS_STAT_MAX;
     }
     psTrace("psLib.math", 4, "---- %s(%d) end ----\n", __func__, numValid);
     return numValid;
 }
-
-#if 0
-/******************************************************************************
-vectorCheckNonEmpty(myVector, maskVector, maskVal, stats): This routine returns
-"true" if the inputPsVector has 1 or more valid elements (a valid element is an
-unmasked element within the specified min/max range).  Otherwise, return
-"false".
- *****************************************************************************/
-static bool vectorCheckNonEmpty(const psVector* myVector,
-                                const psVector* maskVector,
-                                psMaskType maskVal,
-                                psStats* stats)
-{
-    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_VECTOR_NON_NULL(myVector, false);
-    PS_ASSERT_PTR_NON_NULL(stats, false);
-
-    if (stats->options & PS_STAT_USE_RANGE) {
-        if (maskVector) {
-            for (long i = 0; i < myVector->n; i++) {
-                if (!(maskVal & maskVector->data.U8[i]) &&
-                        (stats->min <= myVector->data.F32[i]) &&
-                        (myVector->data.F32[i] <= stats->max)) {
-                    psTrace("psLib.math", 4, "---- %s(true) end ----\n", __func__);
-                    return true;
-                }
-            }
-        } else {
-            for (long i = 0; i < myVector->n; i++) {
-                if ((stats->min <= myVector->data.F32[i]) &&
-                        (myVector->data.F32[i] <= stats->max)) {
-                    psTrace("psLib.math", 4, "---- %s(true) end ----\n", __func__);
-                    return true;
-                }
-            }
-        }
-    } else {
-        if (maskVector != NULL) {
-            for (long i = 0; i < myVector->n; i++) {
-                if (!(maskVal & maskVector->data.U8[i])) {
-                    psTrace("psLib.math", 4, "---- %s(true) end ----\n", __func__);
-                    return true;
-                }
-            }
-        } else {
-            if (myVector->n > 0) {
-                psTrace("psLib.math", 4, "---- %s(true) end ----\n", __func__);
-                return true;
-            }
-        }
-    }
-    psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__);
-    return(false);
-}
-#endif
 
 /******************************************************************************
 vectorSampleMedian(myVector, maskVector, maskVal, stats): calculates the median
 and quartiles of the input vector.  Returns true on success (including if there
-were no valid input vector elements).
- *****************************************************************************/
-static bool vectorSampleMedian(const psVector* myVector,
+were no valid input vector elements). Expects F32 vector for input.
+*****************************************************************************/
+static bool vectorSampleMedian(const psVector* inVector,
                                const psVector* maskVector,
                                psMaskType maskVal,
@@ -428 +240 @@
     psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
 
+    bool useRange = stats->options & PS_STAT_USE_RANGE;
+    psU8 *maskData = (maskVector == NULL) ? NULL : maskVector->data.U8; // Dereference the vector
+    psF32 *input = inVector->data.F32; // Dereference the vector
+
     // Allocate temporary vectors for the data.
-    psVector* vector = psVectorAllocEmpty(myVector->n, PS_TYPE_F32); // Vector containing the unmasked values
+    psVector *outVector = psVectorAllocEmpty(inVector->n, PS_TYPE_F32); // Vector containing the unmasked values
+    psF32 *output = outVector->data.F32; // Dereference the vector
+
     long count = 0;                     // Number of valid entries
 
-    // Determine if we must only use data points within a min/max range.
-    if (stats->options & PS_STAT_USE_RANGE) {
-        // Store all non-masked data points within the min/max range
-        // into the temporary vectors.
-        if (maskVector != NULL) {
-            for (long i = 0; i < myVector->n; i++) {
-                if (!(maskVal & maskVector->data.U8[i]) &&
-                        (stats->min <= myVector->data.F32[i]) &&
-                        (myVector->data.F32[i] <= stats->max)) {
-                    vector->data.F32[count] = myVector->data.F32[i];
-                    count++;
-                }
-            }
-        } else {
-            for (long i = 0; i < myVector->n; i++) {
-                if ((stats->min <= myVector->data.F32[i]) &&
-                        (myVector->data.F32[i] <= stats->max)) {
-                    vector->data.F32[count] = myVector->data.F32[i];
-                    count++;
-                }
-            }
-        }
-    } else {
-        // Store all non-masked data points into the temporary vectors.
-        if (maskVector) {
-            for (long i = 0; i < myVector->n; i++) {
-                if (!(maskVal & maskVector->data.U8[i])) {
-                    vector->data.F32[count] = myVector->data.F32[i];
-                    count++;
-                }
-            }
-        } else {
-            vector = psVectorCopy(vector, myVector, PS_TYPE_F32);
-            count = myVector->n;
-        }
-    }
-    vector->n = count;
+    // Store all non-masked data points within the min/max range
+    // into the temporary vectors.
+    for (long i = 0; i < inVector->n; i++) {
+        if (useRange && (input[i] < stats->min))
+            continue;
+        if (useRange && (input[i] > stats->max))
+            continue;
+        if (maskData && (maskData[i] & maskVal))
+            continue;
+
+        output[count] = input[i];
+        count++;
+    }
+    outVector->n = count;
+
     if (count == 0) {
         psError(PS_ERR_BAD_PARAMETER_SIZE, true, "No valid data in input vector.\n");
-        psFree(vector);
+        psFree(outVector);
+        stats->sampleUQ = NAN;
+        stats->sampleLQ = NAN;
         stats->sampleMedian = NAN;
         return false;
@@ -477 +275 @@
 
     // Sort the temporary vector.
-    vector = psVectorSort(vector, vector); // Sort in-place (since it's a copy, it's OK)
-    if (!vector) {
-        psError(PS_ERR_UNEXPECTED_NULL,
-                false,
-                _("Failed to sort input data."));
-        psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__);
-        psFree(vector);
+    // XXX this is messed up: should psVectorSort free on error or not?
+    if (!psVectorSort(outVector, outVector)) { // Sort in-place (since it's a copy, it's OK)
+        psError(PS_ERR_UNEXPECTED_NULL, false, _("Failed to sort input data."));
+        stats->sampleUQ = NAN;
+        stats->sampleLQ = NAN;
+        stats->sampleMedian = NAN;
         return false;
     }
 
-    // Calculate the median exactly.  Use the average if the number of samples
-    // is even.
+    // Calculate the median.  Use the average if the number of samples if even.
     if (count % 2 == 0) {
-        stats->sampleMedian = 0.5 * (vector->data.F32[(count / 2) - 1] +
-                                     vector->data.F32[count / 2]);
+        stats->sampleMedian = 0.5 * (output[(count/2) - 1] + output[count/2]);
     } else {
-        stats->sampleMedian = vector->data.F32[count / 2];
+        stats->sampleMedian = output[count/2];
     }
 
     // Calculate the quartile points exactly.
-    stats->sampleUQ = vector->data.F32[3 * (count / 4)];
-    stats->sampleLQ = vector->data.F32[count / 4];
+    stats->sampleUQ = output[(int)(0.75*count)];
+    stats->sampleLQ = output[(int)(0.25*count)];
+
+    stats->results |= PS_STAT_SAMPLE_MEDIAN;
+    stats->results |= PS_STAT_SAMPLE_QUARTILE;
 
     // Free the temporary data structures.
-    psFree(vector);
+    psFree(outVector);
 
     // Return "true" on success.
     psTrace("psLib.math", 4, "---- %s(true) end ----\n", __func__);
     return true;
-}
-
-/******************************************************************************
-vectorSmoothHistGaussian(): This routine smoothes the data in the input
-robustHistogram with a Gaussian of width sigma.  It returns a psVector of the
-smoothed data.
- *****************************************************************************/
-psVector *vectorSmoothHistGaussian(psHistogram *histogram,
-                                   psF32 sigma)
-{
-    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
-    psTrace("psLib.math", 5, "(histogram->nums->n, sigma) is (%d, %.2f\n", (int) histogram->nums->n, sigma);
-    PS_ASSERT_PTR_NON_NULL(histogram, NULL);
-    PS_ASSERT_PTR_NON_NULL(histogram->bounds, NULL);
-    PS_ASSERT_PTR_NON_NULL(histogram->nums, NULL);
-    if (psTraceGetLevel("psLib.math") >= 8) {
-        PS_VECTOR_PRINT_F32(histogram->nums);
-    }
-
-    long numBins = histogram->nums->n;  // Number of histogram bins
-    psVector *smooth = psVectorAlloc(numBins, PS_TYPE_F32); // Smoothed version of histogram bins
-    const psVector *bounds = histogram->bounds; // The bounds for the histogram bins
-
-    if (!histogram->uniform) {
-        //
-        // We get here if the histogram is non-uniform.  This code is not tested.
-        // However, it is also not used anywhere, yet.
-        //
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSmoothHistGaussian() on non-uniform histograms has not been tested or used.\n");
-
-        for (long i = 0; i < numBins; i++) {
-            // Determine the midpoint of bin i.
-            float iMid = PS_BIN_MIDPOINT(histogram, i);
-
-            //
-            // We determine the bin numbers (jMin:jMax) corresponding to a
-            // range of data values surrounding iMid.  The range is of size:
-            // 2*PS_GAUSS_WIDTH*sigma
-            long jMin, jMax;
-            psF32 x = iMid - (PS_GAUSS_WIDTH * sigma);
-            for (jMin = i; jMin > 0 && bounds->data.F32[jMin] > x; jMin--)
-                ;
-            x = iMid + (PS_GAUSS_WIDTH * sigma);
-            for (jMax = i; jMax < bounds->n - 1 && bounds->data.F32[jMax + 1] > x; jMax++)
-                ;
-
-            //
-            // Loop from jMin to jMax, computing the gaussian of data i.
-            //
-            smooth->data.F32[i] = 0.0;
-            for (long j = jMin ; j <= jMax ; j++) {
-                float jMid = PS_BIN_MIDPOINT(histogram, j);
-                smooth->data.F32[i] +=
-                    histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true);
-            }
-        }
-    } else {
-        //
-        // We get here if the histogram is uniform.
-        //
-        for (long i = 0; i < numBins; i++) {
-            psF32 binSize = bounds->data.F32[1] - bounds->data.F32[0];
-            psS32 gaussWidth = ((PS_GAUSS_WIDTH * sigma) / binSize);
-
-            //
-            // XXX: The following is wrong.  However, in practice, the sigma was too
-            // large compared to the binSize.  This meant that the smoothing was done
-            // over 500 bins in the robust stats algorithm.  This mean that the smoothing
-            // took way too long.
-            //
-            #if 0
-
-            gaussWidth = 10;
-            #endif
-            //
-            // We determine the bin numbers (jMin:jMax) corresponding to a
-            // range of data values surrounding iMid.  The range is of size:
-            // 2*PS_GAUSS_WIDTH*sigma
-            //
-            psS32 jMin = PS_MAX(i - gaussWidth, 0);
-            psS32 jMax = PS_MIN(i + gaussWidth, bounds->n - 1);
-
-            //
-            // Loop from jMin to jMax, computing the gaussian of data i.
-            //
-            smooth->data.F32[i] = 0.0;
-            float iMid = PS_BIN_MIDPOINT(histogram, i);
-            for (long j = jMin ; j <= jMax ; j++) {
-                float jMid = PS_BIN_MIDPOINT(histogram, j);
-                smooth->data.F32[i] +=
-                    histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true);
-            }
-        }
-    }
-
-    if (psTraceGetLevel("psLib.math") >= 8) {
-        PS_VECTOR_PRINT_F32(smooth);
-    }
-    psTrace("psLib.math", 4, "---- %s() end ----\n", __func__);
-    return(smooth);
 }
 
@@ -619 +317 @@
     NULL
  
- *****************************************************************************/
+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:
+ 
+                     single
+(mask: 0, range: 0): 0.193 sec
+(mask: 1, range: 0): 0.257 sec
+(mask: 0, range: 0): 0.349 sec
+(mask: 1, range: 0): 0.401 sec
+ 
+*****************************************************************************/
 static bool vectorSampleStdev(const psVector* myVector,
                               const psVector* errors,
@@ -630 +337 @@
     // This procedure requires the mean.  If it has not been already
     // calculated, then call vectorSampleMean()
-    if (isnan(stats->sampleMean)) {
+    if (!(stats->results & PS_STAT_SAMPLE_MEAN)) {
         vectorSampleMean(myVector, errors, maskVector, maskVal, stats);
     }
+
     // If the mean is NAN, then generate a warning and set the stdev to NAN.
     if (isnan(stats->sampleMean)) {
@@ -641 +349 @@
     }
 
+    psF32 *data = myVector->data.F32;   // Dereference
+    psU8 *maskData = (maskVector == NULL) ? NULL : maskVector->data.U8;
+    bool useRange = stats->options & PS_STAT_USE_RANGE;
+    psF32 *errorsData = (errors == NULL) ? NULL : errors->data.F32;
+
     // Accumulate the sums
     double mean = stats->sampleMean;    // The mean
@@ -647 +360 @@
     double errorDivisor = 0.0;          // Division for errors
     long count = 0;                     // Number of data points being used
-    if (stats->options & PS_STAT_USE_RANGE) {
-        if (maskVector) {
-            for (long i = 0; i < myVector->n; i++) {
-                if (!(maskVal & maskVector->data.U8[i]) &&
-                        (stats->min <= myVector->data.F32[i]) &&
-                        (myVector->data.F32[i] <= stats->max)) {
-                    double diff = myVector->data.F32[i] - mean;
-                    sumSquares += PS_SQR(diff);
-                    sumDiffs += diff;
-                    count++;
-                    if (errors) {
-                        errorDivisor += 1.0 / PS_SQR(errors->data.F32[i]);
-                    }
-                }
-            }
-        } else {
-            for (long i = 0; i < myVector->n; i++) {
-                if ((stats->min <= myVector->data.F32[i]) && (myVector->data.F32[i] <= stats->max)) {
-                    double diff = myVector->data.F32[i] - mean;
-                    sumSquares += PS_SQR(diff);
-                    sumDiffs += diff;
-                    count++;
-                    if (errors) {
-                        errorDivisor += 1.0 / PS_SQR(errors->data.F32[i]);
-                    }
-                }
-            }
-            count = myVector->n;
-        }
-    } else {
-        if (maskVector) {
-            for (long i = 0; i < myVector->n; i++) {
-                if (!(maskVal & maskVector->data.U8[i])) {
-                    double diff = myVector->data.F32[i] - mean;
-                    sumSquares += PS_SQR(diff);
-                    sumDiffs += diff;
-                    count++;
-                    if (errors) {
-                        errorDivisor += 1.0 / PS_SQR(errors->data.F32[i]);
-                    }
-                }
-            }
-        } else {
-            for (long i = 0; i < myVector->n; i++) {
-                double diff = myVector->data.F32[i] - mean;
-                sumSquares += PS_SQR(diff);
-                sumDiffs += diff;
-                count++;
-                if (errors) {
-                    errorDivisor += 1.0 / PS_SQR(errors->data.F32[i]);
-                }
-            }
-            count = myVector->n;
+    for (long i = 0; i < myVector->n; 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;
+
+        double diff = data[i] - mean;
+        sumSquares += PS_SQR(diff);
+        sumDiffs += diff;
+        count++;
+        if (errors) {
+            errorDivisor += 1.0 / PS_SQR(errorsData[i]);
         }
     }
@@ -717 +392 @@
         stats->sampleStdev = (1.0 / sqrtf(errorDivisor));
     } else {
-        stats->sampleStdev = sqrt((sumSquares - (sumDiffs * sumDiffs / (float)count)) /
-                                  (float)(count - 1));
-    }
+        stats->sampleStdev = sqrt((sumSquares - (sumDiffs * sumDiffs / (float)count)) / (float)(count - 1));
+    }
+    stats->results |= PS_STAT_SAMPLE_STDEV;
+
     psTrace("psLib.math", 4, "---- %s() end ----\n", __func__);
 
@@ -736 +412 @@
 Returns
     true for success; false otherwise
- *****************************************************************************/
+*****************************************************************************/
 static bool vectorClippedStats(const psVector* myVector,
                                const psVector* errors,
@@ -757 +433 @@
                                  PS_CLIPPED_SIGMA_UB, -1);
 
-    // Allocate a psStats structure for calculating the mean, median, and
-    // stdev.
-    psStats *statsTmp = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
+    // unless we succeed, these will have NAN values
+    stats->clippedMean = NAN;
+    stats->clippedStdev = NAN;
+    stats->clippedNvalues = 0;
 
     // We copy the mask vector, to preserve the original
@@ -773 +450 @@
     }
 
-    // 1. Compute the sample median.
-    vectorSampleMedian(myVector, tmpMask, maskVal, statsTmp);
-    if (isnan(statsTmp->sampleMedian)) {
+    // 1. Compute the sample median, which we save for output
+    vectorSampleMedian(myVector, tmpMask, maskVal, stats);
+    if (isnan(stats->sampleMedian)) {
         psLogMsg(__func__, PS_LOG_WARN, "Call to vectorSampleMedian returned NAN\n");
-        stats->clippedMean = NAN;
-        stats->clippedStdev = NAN;
         psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__);
         psFree(tmpMask);
-        psFree(statsTmp);
         return false;
     }
-    psTrace("psLib.math", 6, "The initial sample median is %f\n", statsTmp->sampleMedian);
-
-    // 2. Compute the sample standard deviation.
-    vectorSampleStdev(myVector, errors, tmpMask, maskVal, statsTmp);
-    if (isnan(statsTmp->sampleStdev)) {
+    psTrace("psLib.math", 6, "The initial sample median is %f\n", stats->sampleMedian);
+
+    // 2. Compute the sample standard deviation, which we save for output
+    vectorSampleStdev(myVector, errors, tmpMask, maskVal, stats);
+    if (isnan(stats->sampleStdev)) {
         psLogMsg(__func__, PS_LOG_WARN, "Call to vectorSampleStdev returned NAN\n");
-        stats->clippedMean = NAN;
-        stats->clippedStdev = NAN;
         psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__);
         psFree(tmpMask);
-        psFree(statsTmp);
         return false;
     }
-    psTrace("psLib.math", 6, "The initial sample stdev is %f\n", statsTmp->sampleStdev);
+    psTrace("psLib.math", 6, "The initial sample stdev is %f\n", stats->sampleStdev);
 
     // 3. Use the sample median as the first estimator of the mean X.
-    psF32 clippedMean = statsTmp->sampleMedian;
+    psF32 clippedMean = stats->sampleMedian;
 
     // 4. Use the sample stdev as the first estimator of the mean stdev.
-    psF32 clippedStdev = statsTmp->sampleStdev;
+    psF32 clippedStdev = stats->sampleStdev;
 
     // 5. Repeat N (stats->clipIter) times:
@@ -836 +507 @@
 
         // b) compute new mean and stdev
+        // Allocate a psStats structure for calculating the mean and stdev.
+        psStats *statsTmp = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
         vectorSampleMean(myVector, errors, tmpMask, maskVal, statsTmp);
         vectorSampleStdev(myVector, errors, tmpMask, maskVal, statsTmp);
@@ -853 +526 @@
             clippedStdev = statsTmp->sampleStdev;
         }
+        psFree(statsTmp);
     }
 
@@ -859 +533 @@
 
     // 7. The last calcuated value of x is the clipped mean.
-    if (stats->options & PS_STAT_CLIPPED_MEAN) {
-        stats->clippedMean = clippedMean;
-        psTrace("psLib.math", 6, "The final clipped mean is %f\n", clippedMean);
-    }
     // 8. The last calcuated value of stdev is the clipped stdev.
-    if (stats->options & PS_STAT_CLIPPED_STDEV) {
-        stats->clippedStdev = clippedStdev;
-        psTrace("psLib.math", 6, "The final clipped stdev is %f\n", clippedStdev);
-    }
+    // we always return both stats even if only one was requested
+    stats->clippedMean = clippedMean;
+    stats->clippedStdev = clippedStdev;
+
+    stats->results |= PS_STAT_CLIPPED_MEAN;
+    stats->results |= PS_STAT_CLIPPED_STDEV;
+
+    psTrace("psLib.math", 6, "The final clipped mean is %f\n", clippedMean);
+    psTrace("psLib.math", 6, "The final clipped stdev is %f\n", clippedStdev);
 
     psFree(tmpMask);
-    psFree(statsTmp);
     psTrace("psLib.math", 4, "---- %s(true) end ----\n", __func__);
     return true;
 }
-
-static psF32 QuadraticInverse(psF32 a,
-                              psF32 b,
-                              psF32 c,
-                              psF32 y,
-                              psF32 xLo,
-                              psF32 xHi
-                             )
-{
-    psF64 tmp = sqrt((y - c)/a + (b*b)/(4.0 * a * a));
-
-    psF64 x1 = -b/(2.0*a) + tmp;
-    psF64 x2 = -b/(2.0*a) - tmp;
-
-    #if 0
-
-    psF64 y1 = (a * x1 * x1) + (b * x1) + c;
-    psF64 y2 = (a * x2 * x2) + (b * x2) + c;
-    printf("QuadraticInverse: %fx^2 + %fx + %f\n", a, b, c);
-    printf("QuadraticInverse: y is %f\n", y);
-    printf("QuadraticInverse: (x1, x2) is (%f %f)\n", x1, x2);
-    printf("QuadraticInverse: (y1, y2) is (%f %f)\n", y1, y2);
-    #endif
-
-    if (xLo <= x1 && x1 <= xHi) {
-        return x1;
-    }
-    if (xLo <= x2 && x2 <= xHi) {
-        return x2;
-    }
-    return 0.5 * (xLo + xHi);
-}
-
-
-/******************************************************************************
-fitQuadraticSearchForYThenReturnX(*xVec, *yVec, binNum, yVal): A general
-routine which fits a quadratic to three points and returns the x-value
-corresponding to the input y-value.  This routine takes psVectors of x/y pairs
-as input, and fits a quadratic to the 3 points surrounding element binNum in
-the vectors (the midpoint between element i and i+1 is used for x[i]).  It
-then determines for what value x does that quadratic f(x) = yVal (the input
-parameter).
- 
-*****************************************************************************/
-static psF32 fitQuadraticSearchForYThenReturnX(const psVector *xVec,
-        psVector *yVec,
-        psS32 binNum,
-        psF32 yVal
-                                              )
-{
-    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
-    psTrace("psLib.math", 5, "binNum, yVal is (%d, %f)\n", binNum, yVal);
-    if (psTraceGetLevel("psLib.math") >= 8) {
-        PS_VECTOR_PRINT_F32(xVec);
-        PS_VECTOR_PRINT_F32(yVec);
-    }
-
-    PS_ASSERT_VECTOR_NON_NULL(xVec, NAN);
-    PS_ASSERT_VECTOR_NON_NULL(yVec, NAN);
-    PS_ASSERT_VECTOR_TYPE(xVec, PS_TYPE_F32, NAN);
-    PS_ASSERT_VECTOR_TYPE(yVec, PS_TYPE_F32, NAN);
-    //    PS_ASSERT_VECTORS_SIZE_EQUAL(xVec, yVec, NAN);
-    PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, (int)(xVec->n - 1), NAN);
-    PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, (int)(yVec->n - 1), NAN);
-
-    psVector *x = psVectorAlloc(3, PS_TYPE_F64);
-    psVector *y = psVectorAlloc(3, PS_TYPE_F64);
-    psF32 tmpFloat = 0.0f;
-
-    if ((binNum >= 1) && (binNum < (yVec->n - 2)) && (binNum < (xVec->n - 2))) {
-        // The general case.  We have all three points.
-        x->data.F64[0] = (psF64) (0.5 * (xVec->data.F32[binNum - 1] + xVec->data.F32[binNum]));
-        x->data.F64[1] = (psF64) (0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum+1]));
-        x->data.F64[2] = (psF64) (0.5 * (xVec->data.F32[binNum+1] + xVec->data.F32[binNum+2]));
-        y->data.F64[0] = yVec->data.F32[binNum - 1];
-        y->data.F64[1] = yVec->data.F32[binNum];
-        y->data.F64[2] = yVec->data.F32[binNum + 1];
-        psTrace("psLib.math", 6, "x vec (orig) is (%f %f %f %f)\n", xVec->data.F32[binNum - 1],
-                xVec->data.F32[binNum], xVec->data.F32[binNum+1], xVec->data.F32[binNum+2]);
-        psTrace("psLib.math", 6, "x data is (%f %f %f)\n", x->data.F64[0], x->data.F64[1], x->data.F64[2]);
-        psTrace("psLib.math", 6, "y data is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]);
-
-        //
-        // Ensure that the y value lies within range of the y values.
-        //
-        if (! (((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[2])) ||
-                ((y->data.F64[2] <= yVal) && (yVal <= y->data.F64[0]))) ) {
-            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
-                    _("Specified yVal, %g, is not within y-range, %g to %g."),
-                    (psF64)yVal, y->data.F64[0], y->data.F64[2]);
-        }
-
-        //
-        // Ensure that the y values are monotonic.
-        //
-        if (((y->data.F64[0] < y->data.F64[1]) && !(y->data.F64[1] <= y->data.F64[2])) ||
-                ((y->data.F64[0] > y->data.F64[1]) && !(y->data.F64[1] >= y->data.F64[2]))) {
-            psError(PS_ERR_UNKNOWN, true,
-                    "This routine must be called with monotonically increasing or decreasing data points.\n");
-            psFree(x);
-            psFree(y);
-            psTrace("psLib.math", 5, "---- %s() end ----\n", __func__);
-            return NAN;
-        }
-
-        // Determine the coefficients of the polynomial.
-        psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
-        myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x);
-        if (myPoly == NULL) {
-            psError(PS_ERR_UNEXPECTED_NULL, false,
-                    _("Failed to fit a 1-dimensional polynomial to the three specified data points.  Returning NAN."));
-            psFree(x);
-            psFree(y);
-            psTrace("psLib.math", 5, "---- %s(NAN) end ----\n", __func__);
-            return NAN;
-        }
-        psTrace("psLib.math", 6, "myPoly->coeff[0] is %f\n", myPoly->coeff[0]);
-        psTrace("psLib.math", 6, "myPoly->coeff[1] is %f\n", myPoly->coeff[1]);
-        psTrace("psLib.math", 6, "myPoly->coeff[2] is %f\n", myPoly->coeff[2]);
-        psTrace("psLib.math", 6, "Fitted y vec is (%f %f %f)\n",
-                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[0]),
-                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[1]),
-                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[2]));
-
-        psTrace("psLib.math", 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal);
-        tmpFloat = QuadraticInverse(myPoly->coeff[2], myPoly->coeff[1], myPoly->coeff[0], yVal,
-                                    x->data.F64[0], x->data.F64[2]);
-        psFree(myPoly);
-
-        if (isnan(tmpFloat)) {
-            psError(PS_ERR_UNEXPECTED_NULL,
-                    false, _("Failed to determine the median of the fitted polynomial.  Returning NAN."));
-            psFree(x);
-            psFree(y);
-            psTrace("psLib.math", 5, "---- %s(NAN) end ----\n", __func__);
-            return(NAN);
-        }
-    } else {
-        // These are special cases where the bin is at the beginning or end of the vector.
-        if (binNum == 0) {
-            // We have two points only at the beginning of the vectors x and y.
-            tmpFloat = 0.5 * (xVec->data.F32[binNum] +
-                              xVec->data.F32[binNum + 1]);
-        } else if (binNum == (xVec->n - 1)) {
-            // The special case where we have two points only at the end of
-            // the vectors x and y.
-            // XXX: Is this right?
-            tmpFloat = xVec->data.F32[binNum];
-        } else if (binNum == (xVec->n - 2)) {
-            // XXX: Is this right?
-            tmpFloat = 0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum + 1]);
-        }
-    }
-
-    psTrace("psLib.math", 6, "FIT: return %f\n", tmpFloat);
-    psFree(x);
-    psFree(y);
-
-    psTrace("psLib.math", 5, "---- %s(%f) end ----\n", __func__, tmpFloat);
-    return tmpFloat;
-}
-
-/******************************************************************************
-NOTE: We assume unnormalized gaussians.
- *****************************************************************************/
-static psF32 minimizeLMChi2Gauss1D(psVector *deriv,
-                                   const psVector *params,
-                                   const psVector *coords
-                                  )
-{
-    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_VECTOR_NON_NULL(params, NAN);
-    PS_ASSERT_VECTOR_SIZE(params, (long)2, NAN);
-    PS_ASSERT_VECTOR_TYPE(params, PS_TYPE_F32, NAN);
-    PS_ASSERT_VECTOR_NON_NULL(coords, NAN);
-    PS_ASSERT_VECTOR_SIZE(coords, (long)1, NAN);
-    PS_ASSERT_VECTOR_TYPE(coords, PS_TYPE_F32, NAN);
-
-    psF32 x = coords->data.F32[0];
-    psF32 mean = params->data.F32[0];
-    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 = dx * gauss;
-        deriv->data.F32[0] = tmp / var;
-        deriv->data.F32[1] = tmp * dx / (var * var);
-    }
-
-
-    psTrace("psLib.math", 4, "---- %s() end ----\n", __func__);
-    return gauss;
-}
-
-/*
- * 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;
-
-    psF32 binSize = 1;
-    if (stats->options & PS_STAT_USE_BINSIZE) {
-        // Set initial bin size to the specified value.
-        binSize = stats->binsize;
-        psTrace("psLib.math", 6, "Setting initial robust bin size to %.2f\n", binSize);
-    } else {
-        // 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));
-        binSize = 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);
-        psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
-        return false;
-    }
-
-    // Calculate the number of bins.
-    long numBins = (max - min) / binSize;
-    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", binSize);
-    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);
-        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(histogram);
-        psFree(x);
-        psFree(y);
-        psFree(minimizer);
-        psFree(params);
-        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;
-}
-
 
 /******************************************************************************
@@ -1278 +596 @@
     }
 
+    // statsMinMax is only applied to a subset of the data points
     psStats *statsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); // Statistics for min and max
     psHistogram *histogram = NULL;      // Histogram of the data
@@ -1287 +606 @@
     // Iterate to get the best bin size
     for (int iterate = 1; iterate > 0; iterate++) {
-        psTrace("psLib.math", 6,
-                "-------------------- Iterating on Bin size.  Iteration number %d --------------------\n",
-                iterate);
+        psTrace("psLib.math", 6, "-------------------- Iterating on Bin size.  Iteration number %d --------------------\n", iterate);
 
         // Get the minimum and maximum values
@@ -1308 +625 @@
         if (fabs(max - min) <= FLT_EPSILON) {
             psTrace("psLib.math", 5, "All data points have the same value: %f.\n", min);
-            if (stats->options & PS_STAT_ROBUST_MEDIAN) {
-                stats->robustMedian = min;
-            }
-            if (stats->options & PS_STAT_ROBUST_QUARTILE) {
-                stats->robustUQ = min;
-                stats->robustLQ = min;
-            }
+            stats->robustMedian = min;
+            stats->robustUQ = min;
+            stats->robustLQ = min;
             stats->robustN50 = numValid;
             psFree(statsMinMax);
@@ -1433 +746 @@
             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
-        // XXX: I am overriding the ADD for now.  The following code implements the ADD exactly and is having
-        // problems fitting a 2nd-order polynomial to data y-values suchs as (1, 1, 100).  Therefore, I
-        // commented out the code and am doing simply linear interpolation.
-
-        // Value for the 15.8655% mark
-        float binLoF32 = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums, binLo,
-                         totalDataPoints * 0.158655f);
-        // Value for the 84.1345% mark
-        float binHiF32 = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums, binHi,
-                         totalDataPoints * 0.841345f);
-        psTrace("psLib.math", 6, "The exact 15.8655%% and 84.1345%% data point positions are: (%f, %f)\n",
-                binLoF32, binHiF32);
-        #else
-        // This code basically interpolates to find the positions exactly.
+
+        // ADD step 4b: Interpolate Sigma (linearly) to find these two positions exactly: these are the 1sigma positions.
         psTrace("psLib.math", 6, "binLo is %ld.  Nums at that bin and the next are (%.2f, %.2f)\n",
                 binLo, cumulative->nums->data.F32[binLo], cumulative->nums->data.F32[binLo+1]);
@@ -1455 +753 @@
                 binHi, cumulative->nums->data.F32[binHi], cumulative->nums->data.F32[binHi+1]);
 
-        float deltaBounds = cumulative->bounds->data.F32[binLo+1] - cumulative->bounds->data.F32[binLo];
-        float deltaNums;
-        float prevPixels;
+        float deltaBounds, deltaNums, prevPixels;
+        float percentNums, base, binLoF32;
+
+        // find the -1 sigma points with linear interpolation
+        deltaBounds = cumulative->bounds->data.F32[binLo+1] - cumulative->bounds->data.F32[binLo];
         if (binLo == 0) {
             deltaNums = cumulative->nums->data.F32[0];
@@ -1465 +765 @@
             prevPixels = cumulative->nums->data.F32[binLo-1];
         }
-        float percentNums = (totalDataPoints * 0.158655f) - prevPixels;
-        float base = cumulative->bounds->data.F32[binLo];
-        float binLoF32 = base + (deltaBounds / deltaNums) * percentNums; // Value for the 15.8655% mark
+        percentNums = (totalDataPoints * 0.158655f) - prevPixels;
+        base = cumulative->bounds->data.F32[binLo];
+        binLoF32 = base + (deltaBounds / deltaNums) * percentNums; // Value for the 15.8655% mark
         psTrace("psLib.math", 6,
                 "(base, deltaBounds, deltaNums, prevPixels, percentNums) is (%.2f %.2f %.2f %.2f %.2f)\n",
                 base, deltaBounds, deltaNums, prevPixels, percentNums);
 
+        // find the +1 sigma points with linear interpolation
         deltaBounds = cumulative->bounds->data.F32[binHi+1] - cumulative->bounds->data.F32[binHi];
         if (binHi == 0) {
@@ -1486 +787 @@
                 "(base, deltaBounds, deltaNums, prevPixels, percentNums) is (%.2f %.2f %.2f %.2f %.2f)\n",
                 base, deltaBounds, deltaNums, prevPixels, percentNums);
-        psTrace("psLib.math", 6,
+
+        // report +/- 1 sigma points
+        psTrace("psLib.math", 5,
                 "The exact 15.8655 and 84.1345 percent data point positions are: (%f, %f)\n",
                 binLoF32, binHiF32);
-        #endif
 
         // ADD step 5: Determine SIGMA as 1/2 of the distance between these positions.
@@ -1585 +887 @@
     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) {
-        stats->fittedStdev = stats->robustStdev;  // pass the guess sigma
-        stats->fittedMean = stats->robustMedian;  // pass the guess mean
-        if (!vectorFittedStats(myVector, errors, mask, maskVal, stats)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to estimate statistics (pass 1)");
-            psFree(statsMinMax);
-            psFree(mask);
-
-            return false;
-        }
-        // if there is a large swing in sigma, try a second time
-        if ((stats->fittedStdev/stats->robustStdev) < 0.75) {
-            if (!vectorFittedStats (myVector, errors, mask, maskVal, stats)) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to estimate statistics (pass 2)");
-                psFree(statsMinMax);
-                psFree(mask);
-
-                return false;
-            }
-        }
-    }
-
     // Clean up
     psFree(statsMinMax);
     psFree(mask);
 
+    stats->results |= PS_STAT_ROBUST_MEDIAN;
+    stats->results |= PS_STAT_ROBUST_STDEV;
+
     psTrace("psLib.math", 4, "---- %s(0) end  ----\n", __func__);
     return true;
 }
 
+/*
+ * 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)
+{
+
+    // This procedure requires the mean.  If it has not been already
+    // calculated, then call vectorSampleMean()
+    if (!(stats->results & PS_STAT_ROBUST_MEDIAN)) {
+        vectorRobustStats(myVector, errors, mask, maskVal, stats);
+    }
+
+    // If the mean is NAN, then generate a warning and set the stdev to NAN.
+    if (isnan(stats->robustMedian)) {
+        stats->fittedStdev = NAN;
+        stats->fittedStdev = NAN;
+        psTrace("psLib.math", 4, "---- %s() end ----\n", __func__);
+        return false;
+    }
+
+    float guessStdev = stats->robustStdev;  // pass the guess sigma
+    float guessMean = stats->robustMedian;  // pass the guess mean
+
+    psTrace("psLib.math", 6, "The guess mean  is %f.\n", guessMean);
+    psTrace("psLib.math", 6, "The guess stdev is %f.\n", guessStdev);
+
+    bool done = false;
+    for (int iteration = 0; !done && (iteration < 2); iteration ++) {
+        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;
+
+        psF32 binSize = 1;
+        if (stats->options & PS_STAT_USE_BINSIZE) {
+            // Set initial bin size to the specified value.
+            binSize = stats->binsize;
+            psTrace("psLib.math", 6, "Setting initial robust bin size to %.2f\n", binSize);
+        } else {
+            // 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));
+            binSize = guessStdev / 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);
+            psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
+            return false;
+        }
+
+        // Calculate the number of bins.
+        // XXX can we calculate the binMin, binMax **before** building this histogram?
+        long numBins = (max - min) / binSize;
+        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", binSize);
+        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 = guessMean - minFitSigma*guessStdev;
+        if (tmpScalar.data.F32 < min) {
+            binMin = 0;
+        } else {
+            binMin = p_psVectorBinDisect(histogram->bounds, &tmpScalar);
+        }
+        tmpScalar.data.F32 = guessMean + maxFitSigma*guessStdev;
+        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);
+            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);
+        psTrace("psLib.math", 6, "The clipped min is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMin), binMin);
+        psTrace("psLib.math", 6, "The clipped max is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMax), binMax);
+
+        // 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] = guessMean;
+        params->data.F32[1] = PS_SQR(guessStdev);
+        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(params);
+            psFree(minimizer);
+            psFree(x);
+            psFree(y);
+            psFree(histogram);
+            psFree(statsMinMax);
+            psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
+            return false;
+        }
+        if (psTraceGetLevel("psLib.math") >= 8) {
+            PS_VECTOR_PRINT_F32(params);
+        }
+
+        guessMean = params->data.F32[0];
+        guessStdev = sqrt(params->data.F32[1]);
+        if (guessStdev > 0.75*stats->robustStdev) {
+            done = true;
+        }
+
+        // Clean up after fitting
+        psFree (params);
+        psFree (minimizer);
+        psFree (x);
+        psFree (y);
+        psFree (histogram);
+        psFree (statsMinMax);
+    }
+
+    // The fitted mean is the Gaussian mean.
+    stats->fittedMean = guessMean;
+    psTrace("psLib.math", 6, "The fitted mean is %f.\n", stats->fittedMean);
+
+    // The fitted standard deviation
+    stats->fittedStdev = guessStdev;
+    psTrace("psLib.math", 6, "The fitted stdev is %f.\n", stats->fittedStdev);
+
+    stats->results |= PS_STAT_FITTED_MEAN;
+    stats->results |= PS_STAT_FITTED_STDEV;
+
+    return true;
+}
+
+
+/********************
+ * vectorFittedStats_v2 requires guess for fittedMean and fittedStdev
+ * robustN50 should also be set
+ * gaussian fit is performed using 2D polynomial to ln(y)
+ ********************/
+static bool vectorFittedStats_v2 (const psVector* myVector,
+                                  const psVector* errors,
+                                  psVector* mask,
+                                  psMaskType maskVal,
+                                  psStats* stats)
+{
+
+    // This procedure requires the mean.  If it has not been already
+    // calculated, then call vectorSampleMean()
+    if (!(stats->results & PS_STAT_ROBUST_MEDIAN)) {
+        vectorRobustStats(myVector, errors, mask, maskVal, stats);
+    }
+
+    // If the mean is NAN, then generate a warning and set the stdev to NAN.
+    if (isnan(stats->robustMedian)) {
+        stats->fittedStdev = NAN;
+        stats->fittedStdev = NAN;
+        psTrace("psLib.math", 4, "---- %s() end ----\n", __func__);
+        return false;
+    }
+
+    float guessStdev = stats->robustStdev;  // pass the guess sigma
+    float guessMean = stats->robustMedian;  // pass the guess mean
+
+    psTrace("psLib.math", 6, "The guess mean  is %f.\n", guessMean);
+    psTrace("psLib.math", 6, "The guess stdev is %f.\n", guessStdev);
+
+    bool done = false;
+    for (int iteration = 0; !done && (iteration < 2); iteration ++) {
+        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;
+
+        psF32 binSize = 1;
+        if (stats->options & PS_STAT_USE_BINSIZE) {
+            // Set initial bin size to the specified value.
+            binSize = stats->binsize;
+            psTrace("psLib.math", 6, "Setting initial robust bin size to %.2f\n", binSize);
+        } else {
+            // 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));
+            binSize = guessStdev / 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);
+            psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
+            return false;
+        }
+
+        // Calculate the number of bins.
+        // XXX can we calculate the binMin, binMax **before** building this histogram?
+        long numBins = (max - min) / binSize;
+        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", binSize);
+        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 = guessMean - minFitSigma*guessStdev;
+        if (tmpScalar.data.F32 < min) {
+            binMin = 0;
+        } else {
+            binMin = p_psVectorBinDisect(histogram->bounds, &tmpScalar);
+        }
+        tmpScalar.data.F32 = guessMean + maxFitSigma*guessStdev;
+        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);
+            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 = psVectorAllocEmpty((1 + (binMax - binMin)), PS_TYPE_F32); // Vector of coordinates
+        psVector *x = psVectorAllocEmpty((1 + (binMax - binMin)), PS_TYPE_F32); // Vector of ordinates
+        long j = 0;
+        for (long i = binMin; i <= binMax ; i++) {
+            if (histogram->nums->data.F32[i] <= 0.0)
+                continue;
+            x->data.F32[j] = PS_BIN_MIDPOINT(histogram, i);
+            y->data.F32[j] = log(histogram->nums->data.F32[i]);
+            // XXX note this is the natural log: expected distribution is A exp(-(x-xo)^2/2sigma^2)
+            j++;
+        }
+        y->n = x->n = j;
+        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);
+        psTrace("psLib.math", 6, "The clipped min is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMin), binMin);
+        psTrace("psLib.math", 6, "The clipped max is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMax), binMax);
+
+        // fit 2nd order polynomial to ln(y) = -(x-xo)^2/2sigma^2
+        psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
+
+        // XXX how can we protect against some extreme outliers?  the robust histogram
+        // should have already dealt with those, but we are again sensitive to them...
+        // psStats *fitStats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+        // fitStats->clipIter = 3.0;
+        // fitStats->clipSigma = 3.0;
+        // psVector *fitMask = psVectorAlloc(y->n, PS_TYPE_U8);
+        // psVectorInit (fitMask, 0);
+
+        if (!psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x)) {
+            psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n");
+            psFree(x);
+            psFree(y);
+            psFree(poly);
+            // psFree(fitStats);
+            // psFree(fitMask);
+            psFree(histogram);
+            psFree(statsMinMax);
+            psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
+            return false;
+        }
+
+        guessStdev = sqrt(-0.5/poly->coeff[2]);
+        guessMean = poly->coeff[1]*PS_SQR(guessStdev);
+        if (guessStdev > 0.75*stats->robustStdev) {
+            done = true;
+        } else {
+            psTrace("psLib.math", 6, "The new guess mean  is %f.\n", guessMean);
+            psTrace("psLib.math", 6, "The new guess stdev is %f.\n", guessStdev);
+        }
+
+        // Clean up after fitting
+        psFree (x);
+        psFree (y);
+        psFree (poly);
+        // psFree (fitStats);
+        // psFree (fitMask);
+        psFree (histogram);
+        psFree (statsMinMax);
+    }
+
+    // The fitted mean is the Gaussian mean.
+    stats->fittedMean = guessMean;
+    psTrace("psLib.math", 6, "The fitted mean is %f.\n", stats->fittedMean);
+
+    // The fitted standard deviation
+    stats->fittedStdev = guessStdev;
+    psTrace("psLib.math", 6, "The fitted stdev is %f.\n", stats->fittedStdev);
+
+    stats->results |= PS_STAT_FITTED_MEAN_V2;
+    stats->results |= PS_STAT_FITTED_STDEV_V2;
+
+    return true;
+}
+
+
+/******************************************************************************
+vectorSmoothHistGaussian(): This routine smoothes the data in the input
+robustHistogram with a Gaussian of width sigma.  It returns a psVector of the
+smoothed data.
+ 
+XXX this function is unused
+*****************************************************************************/
+psVector *vectorSmoothHistGaussian(psHistogram *histogram,
+                                   psF32 sigma)
+{
+    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
+    psTrace("psLib.math", 5, "(histogram->nums->n, sigma) is (%d, %.2f\n", (int) histogram->nums->n, sigma);
+    PS_ASSERT_PTR_NON_NULL(histogram, NULL);
+    PS_ASSERT_PTR_NON_NULL(histogram->bounds, NULL);
+    PS_ASSERT_PTR_NON_NULL(histogram->nums, NULL);
+    if (psTraceGetLevel("psLib.math") >= 8) {
+        PS_VECTOR_PRINT_F32(histogram->nums);
+    }
+
+    long numBins = histogram->nums->n;  // Number of histogram bins
+    psVector *smooth = psVectorAlloc(numBins, PS_TYPE_F32); // Smoothed version of histogram bins
+    const psVector *bounds = histogram->bounds; // The bounds for the histogram bins
+
+    if (!histogram->uniform) {
+        //
+        // We get here if the histogram is non-uniform.  This code is not tested.
+        // However, it is also not used anywhere, yet.
+        //
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSmoothHistGaussian() on non-uniform histograms has not been tested or used.\n");
+
+        for (long i = 0; i < numBins; i++) {
+            // Determine the midpoint of bin i.
+            float iMid = PS_BIN_MIDPOINT(histogram, i);
+
+            //
+            // We determine the bin numbers (jMin:jMax) corresponding to a
+            // range of data values surrounding iMid.  The range is of size:
+            // 2*PS_GAUSS_WIDTH*sigma
+            long jMin, jMax;
+            psF32 x = iMid - (PS_GAUSS_WIDTH * sigma);
+            for (jMin = i; jMin > 0 && bounds->data.F32[jMin] > x; jMin--)
+                ;
+            x = iMid + (PS_GAUSS_WIDTH * sigma);
+            for (jMax = i; jMax < bounds->n - 1 && bounds->data.F32[jMax + 1] > x; jMax++)
+                ;
+
+            //
+            // Loop from jMin to jMax, computing the gaussian of data i.
+            //
+            smooth->data.F32[i] = 0.0;
+            for (long j = jMin ; j <= jMax ; j++) {
+                float jMid = PS_BIN_MIDPOINT(histogram, j);
+                smooth->data.F32[i] +=
+                    histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true);
+            }
+        }
+    } else {
+        //
+        // We get here if the histogram is uniform.
+        //
+        for (long i = 0; i < numBins; i++) {
+            psF32 binSize = bounds->data.F32[1] - bounds->data.F32[0];
+            psS32 gaussWidth = ((PS_GAUSS_WIDTH * sigma) / binSize);
+
+            //
+            // XXX: The following is wrong.  However, in practice, the sigma was too
+            // large compared to the binSize.  This meant that the smoothing was done
+            // over 500 bins in the robust stats algorithm.  This mean that the smoothing
+            // took way too long.
+            //
+            #if 0
+
+            gaussWidth = 10;
+            #endif
+            //
+            // We determine the bin numbers (jMin:jMax) corresponding to a
+            // range of data values surrounding iMid.  The range is of size:
+            // 2*PS_GAUSS_WIDTH*sigma
+            //
+            psS32 jMin = PS_MAX(i - gaussWidth, 0);
+            psS32 jMax = PS_MIN(i + gaussWidth, bounds->n - 1);
+
+            //
+            // Loop from jMin to jMax, computing the gaussian of data i.
+            //
+            smooth->data.F32[i] = 0.0;
+            float iMid = PS_BIN_MIDPOINT(histogram, i);
+            for (long j = jMin ; j <= jMax ; j++) {
+                float jMid = PS_BIN_MIDPOINT(histogram, j);
+                smooth->data.F32[i] +=
+                    histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true);
+            }
+        }
+    }
+
+    if (psTraceGetLevel("psLib.math") >= 8) {
+        PS_VECTOR_PRINT_F32(smooth);
+    }
+    psTrace("psLib.math", 4, "---- %s() end ----\n", __func__);
+    return(smooth);
+}
+
 /*****************************************************************************/
 
@@ -1621 +1395 @@
 
 /*****************************************************************************/
-
-static void histogramFree(psHistogram* myHist);
 
 // We keep statsFree so that we can identify statistics pointers from the memblock
@@ -1632 +1404 @@
 /******************************************************************************
     psStatsAlloc(): This routine must create a new psStats data structure.
- *****************************************************************************/
+*****************************************************************************/
 psStats* psStatsAlloc(psStatsOptions options)
 {
@@ -1675 +1447 @@
 }
 
-
-
-/******************************************************************************
-psHistogramAlloc(lower, upper, n): allocate a uniform histogram structure
-with the specifed upper and lower limits, and the specifed number of bins.
-This routine will also set the bounds for each of the bins.
- 
-Input:
-    lower
-    upper
-    n
-Returns:
-    The histogram structure
- *****************************************************************************/
-psHistogram* psHistogramAlloc(float lower, float upper, int n)
-{
-    psTrace("psLib.math", 3, "---- %s() begin  ----\n", __func__);
-    psTrace("psLib.math", 5, "(lower, upper, n) is (%f, %f, %d)\n", lower, upper, n);
-    PS_ASSERT_INT_POSITIVE(n, NULL);
-    PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(upper, lower, NULL);
-
-    // Allocate memory for the new histogram structure.  If there are N bins, then there are N+1 bounds to
-    // those bins.
-    psHistogram *newHist = (psHistogram* ) psAlloc(sizeof(psHistogram)); // The new histogram structure
-    psMemSetDeallocator(newHist, (psFreeFunc) histogramFree);
-    psVector* newBounds = psVectorAlloc(n + 1, PS_TYPE_F32);
-    newHist->bounds = newBounds;
-
-    // Calculate the bounds for each bin.
-    psF32 binSize = (upper - lower) / (psF32)n; // The histogram bin size
-    // XXX: Is the following necessary? It prevents the max data point from being in a non-existant bin.
-    binSize += FLT_EPSILON;
-    for (long i = 0; i < n + 1; i++) {
-        newBounds->data.F32[i] = lower + (binSize * (psF32)i);
-    }
-
-    // Allocate the bins, and initialize them to zero.
-    newHist->nums = psVectorAlloc(n, PS_TYPE_F32);
-    psVectorInit(newHist->nums, 0.0);
-
-    // Initialize the other members.
-    newHist->minNum = 0;
-    newHist->maxNum = 0;
-    newHist->uniform = true;
-
-    psTrace("psLib.math", 3, "---- %s() end  ----\n", __func__);
-    return newHist;
-}
-
-/******************************************************************************
-psHistogramAllocGeneric(bounds): allocate a non-uniform histogram structure
-with the specifed bounds.
- 
-Input:
-    bounds
-Returns:
-    The histogram structure
- *****************************************************************************/
-psHistogram* psHistogramAllocGeneric(const psVector* bounds)
-{
-    psTrace("psLib.math", 3, "---- %s() begin  ----\n", __func__);
-    PS_ASSERT_VECTOR_NON_NULL(bounds, NULL);
-    PS_ASSERT_VECTOR_TYPE(bounds, PS_TYPE_F32, NULL);
-    PS_ASSERT_LONG_LARGER_THAN_OR_EQUAL(bounds->n, (long)2, NULL);
-
-    // Allocate memory for the new histogram structure.
-    psHistogram *newHist = (psHistogram* ) psAlloc(sizeof(psHistogram)); // The new histogram structure
-    psMemSetDeallocator(newHist, (psFreeFunc) histogramFree);
-    psVector* newBounds = psVectorCopy(NULL, bounds, PS_TYPE_F32);
-    newHist->bounds = newBounds;
-
-    // Allocate the bins, and initialize them to zero.  If there are N bounds,
-    // then there are N-1 bins.
-    newHist->nums = psVectorAlloc((bounds->n) - 1, PS_TYPE_F32);
-    psVectorInit(newHist->nums, 0.0);
-
-    // Initialize the other members.
-    newHist->minNum = 0;
-    newHist->maxNum = 0;
-    newHist->uniform = false;
-
-    psTrace("psLib.math", 3, "---- %s() end  ----\n", __func__);
-    return (newHist);
-}
-
-static void histogramFree(psHistogram* myHist)
-{
-    psFree(myHist->bounds);
-    psFree(myHist->nums);
-}
-
-
-bool psMemCheckHistogram(psPtr ptr)
-{
-    PS_ASSERT_PTR(ptr, false);
-    return ( psMemGetDeallocator(ptr) == (psFreeFunc)histogramFree );
-}
-
-
-
-/*****************************************************************************
-UpdateHistogramBins(binNum, out, data, error): This routine is to be used when
-updating the histogram in the presence of errors in the input data.  We treat
-the data point as a boxcar PDF and update a range of points surrounding the
-histogram bin which contains the point.  The width of that boxcar is defined
-as 2.35 * error.
- 
-XXX: Must test this.
- *****************************************************************************/
-static bool UpdateHistogramBins(long binNum, // Bin number of the data point
-                                psHistogram* out, // The histogram to be updated
-                                psF32 data, // The data point value
-                                psF32 error // The error in the data point
-                               )
-{
-    psTrace("psLib.math", 3, "---- %s() begin  ----\n", __func__);
-    PS_ASSERT_PTR_NON_NULL(out, false);
-    PS_ASSERT_PTR_NON_NULL(out->bounds, false);
-    PS_ASSERT_PTR_NON_NULL(out->nums, false);
-    PS_ASSERT_LONG_WITHIN_RANGE(binNum, (long)0, (long)((out->nums->n)-1), false);
-    PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(error, 0.0, false);
-    PS_ASSERT_FLOAT_WITHIN_RANGE(data, out->bounds->data.F32[0],
-                                 out->bounds->data.F32[(out->bounds->n)-1], false);
-
-    psF32 boxcarWidth = 2.35 * error;   // Width of the boxcar
-    psF32 boxcarCenter = (out->bounds->data.F32[binNum] +
-                          out->bounds->data.F32[binNum+1]) / 2.0; // Centre of the boxcar
-    psF32 boxcarLeft = boxcarCenter - (boxcarWidth / 2.0); // Left endpoint of the boxcar for the PDF
-    psF32 boxcarRight = boxcarCenter + (boxcarWidth / 2.0); // Right endpoint of the boxcar for the PDF
-    psS32 boxcarLeftBinNum = 0;         // Bin number for left endpoint
-    psS32 boxcarRightBinNum = 0;        // Bin number for right endpoint
-
-    // Determine the left endpoint of the boxcar for the PDF.
-    for (long bin = binNum; bin >= 0; bin--) {
-        if (out->nums->data.F32[bin] <= boxcarLeft) {
-            boxcarLeftBinNum = bin;
-            break;
-        }
-    }
-
-    // Determine the right endpoint of the boxcar for the PDF.
-    for (long bin = binNum; bin < out->nums->n; bin++) {
-        if (out->nums->data.F32[bin] >= boxcarRight) {
-            boxcarRightBinNum = bin;
-            break;
-        }
-    }
-
-    // If the boxcar fits entirely inside this bin, then simply add 1.0 to the
-    // bin and return.
-    if (boxcarLeftBinNum == boxcarRightBinNum) {
-        out->nums->data.F32[binNum] += 1.0;
-        psTrace("psLib.math", 3, "---- %s(true) end  ----\n", __func__);
-        return true;
-    }
-
-    // If we get here, multiple bins must be updated.  We handle the left-most endpoint, and right-most
-    // endpoints differently.
-    out->nums->data.F32[boxcarLeftBinNum] +=
-        (out->bounds->data.F32[boxcarLeftBinNum+1] - boxcarLeft) / boxcarWidth;
-
-    // Loop through the center bins, if any.
-    for (long bin = boxcarLeftBinNum + 1; bin < (boxcarRightBinNum - 1); bin++) {
-        out->nums->data.F32[bin] +=
-            (out->bounds->data.F32[bin+1] - out->bounds->data.F32[bin]) / boxcarWidth;
-    }
-
-    // Handle the right endpoint differently.
-    out->nums->data.F32[boxcarRightBinNum]+=
-        (boxcarRight - out->bounds->data.F32[boxcarRightBinNum]) / boxcarWidth;
-
-    psTrace("psLib.math", 3, "---- %s(true) end  ----\n", __func__);
-    return true;
-}
-
-
-/*****************************************************************************
-psVectorHistogram(out, in, errors, mask, maskVal): this procedure takes as
-input a preallocated and initialized histogram structure.  It fills the bins
-in that histogram structure in accordance with the input data "in" and the,
-possibly NULL, mask vector.
- 
-Inputs:
-    out
-    in
-    mask
-    maskVal
-Returns:
-    The histogram structure "out".
- *****************************************************************************/
-psHistogram* psVectorHistogram(psHistogram* out,
-                               const psVector* values,
-                               const psVector* errors,
-                               const psVector* mask,
-                               psMaskType maskVal)
-{
-    psTrace("psLib.math", 3, "---- %s() begin  ----\n", __func__);
-    PS_ASSERT_PTR_NON_NULL(out, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(out->bounds, NULL);
-    PS_ASSERT_VECTOR_TYPE(out->bounds, PS_TYPE_F32, NULL);
-    PS_ASSERT_INT_NONNEGATIVE(out->bounds->n, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(out->nums, NULL);
-    PS_ASSERT_VECTOR_TYPE(out->nums, PS_TYPE_F32, NULL);
-    PS_ASSERT_INT_NONNEGATIVE(out->nums->n, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(values, out);
-    if (mask) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(values, mask, NULL);
-        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
-    }
-    if (errors) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(values, errors, NULL);
-        PS_ASSERT_VECTOR_TYPE(errors, values->type.type, NULL);
-    }
-
-    long binNum = 0;                    // A temporary bin number
-    long numBins = out->nums->n;        // The total number of bins
-    psScalar tmpScalar;
-    tmpScalar.type.type = PS_TYPE_F32;
-
-    // Convert input and errors vectors to F32 if necessary.
-    psVector* inF32 = NULL;             // F32 version of input vector
-    if (values->type.type == PS_TYPE_F32) {
-        inF32 = psMemIncrRefCounter((psPtr)values);
-    } else {
-        inF32 = psVectorCopy(NULL, values, PS_TYPE_F32);
-    }
-    psVector* errorsF32 = NULL;         // F32 version of errors vector
-    if (errors) {
-        if (errors->type.type == PS_TYPE_F32) {
-            errorsF32 = psMemIncrRefCounter((psPtr)errors);
-        } else {
-            errorsF32 = psVectorCopy(NULL, errors, PS_TYPE_F32);
-        }
-    }
-
-    for (long i = 0; i < inF32->n; i++) {
-        // Check if this pixel is masked, and if so, skip it.
-        if (!mask || (mask && (!(mask->data.U8[i] & maskVal)))) {
-            if (inF32->data.F32[i] < out->bounds->data.F32[0]) {
-                // If this pixel is below minimum value, count it, then skip.
-                out->minNum++;
-            } else if (inF32->data.F32[i] > out->bounds->data.F32[numBins]) {
-                // If this pixel is above maximum value, count it, then skip.
-                out->maxNum++;
-            } else {
-                // If this is a uniform histogram, determining the correct
-                // number is trivial.
-                if (out->uniform == true) {
-                    float binSize = out->bounds->data.F32[1] - out->bounds->data.F32[0]; // Histogram bin size
-                    binNum = (inF32->data.F32[i] - out->bounds->data.F32[0]) / binSize;
-                    if (errorsF32) {
-                        if (!UpdateHistogramBins(binNum, out, inF32->data.F32[i], errorsF32->data.F32[i])) {
-                            psLogMsg(__func__, PS_LOG_WARN, "WARNING: Failed to update the histogram "
-                                     "bins with the errors vector.\n");
-                        }
-                    } else {
-                        // This if-statement really shouldn't be necessary.
-                        // However, due to numerical lack of precision, we
-                        // occasionally produce a binNum outside the range.
-                        if (binNum >= out->nums->n) {
-                            binNum = out->nums->n - 1;
-                        }
-                        (out->nums->data.F32[binNum])+= 1.0;
-                    }
-
-                } else {
-                    // If this is a non-uniform histogram, determining the
-                    // correct bin number requires a bit more work.
-                    tmpScalar.data.F32 = inF32->data.F32[i];
-                    binNum = p_psVectorBinDisect( *(psVector* *)&out->bounds, &tmpScalar);
-                    if (binNum < 0) {
-                        psLogMsg(__func__, PS_LOG_WARN,
-                                 "WARNING: psVectorHistogram(): element outside histogram bounds.\n");
-                    } else {
-                        if (errorsF32 != NULL) {
-                            if (!UpdateHistogramBins(binNum, out, inF32->data.F32[i], errors->data.F32[i])) {
-                                psLogMsg(__func__, PS_LOG_WARN, "WARNING: Failed to update the histogram "
-                                         "bins with the errors vector.\n");
-                            }
-                        } else {
-                            out->nums->data.F32[binNum] += 1.0;
-                        }
-                    }
-                }
-            }
-        }
-    }
-
-    psFree(inF32);
-    psFree(errorsF32);
-
-    psTrace("psLib.math", 3, "---- %s() end  ----\n", __func__);
-    return (out);
-}
-
 /******************************************************************************
 psVectorStats(in, mask, maskVal, stats): this is the public API
@@ -1984 +1461 @@
  
 XXX: Should we free stats if the asserts fail? NO; we don't own it (RHL).
- *****************************************************************************/
-psStats* psVectorStats(
-    psStats* stats,
-    const psVector* in,
-    const psVector* errors,
-    const psVector* mask,
-    psMaskType maskVal)
+*****************************************************************************/
+bool psVectorStats(psStats* stats,
+                   const psVector* in,
+                   const psVector* errors,
+                   const psVector* mask,
+                   psMaskType maskVal)
 {
     psTrace("psLib.math", 3,"---- %s() begin  ----\n", __func__);
@@ -2036 +1512 @@
     }
 
+    // init the value of stats->results: this is used internally to check if
+    // prior functions have been called
+    stats->results = PS_STAT_NONE;
+    bool status = true;
+
     // ************************************************************************
     if (stats->options & PS_STAT_SAMPLE_MEAN) {
         if (!vectorSampleMean(inF32, errorsF32, maskU8, maskVal, stats)) {
-            psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSampleMean() returned an error.\n");
-            stats->sampleMean = NAN;
+            psError(PS_ERR_UNKNOWN, false, "Failed to calculate vector sample mean");
+            status &= false;
         }
     }
@@ -2047 +1528 @@
     if (stats->options & (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_QUARTILE)) {
         if (!vectorSampleMedian(inF32, maskU8, maskVal, stats)) {
-            psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSampleMedian() returned an error.\n");
-            stats->sampleMedian = NAN;
-            stats->sampleUQ = NAN;
-            stats->sampleLQ = NAN;
+            psError(PS_ERR_UNKNOWN, false, "Failed to calculate sample median");
+            status &= false;
         }
     }
@@ -2056 +1535 @@
     // ************************************************************************
     if (stats->options & PS_STAT_SAMPLE_STDEV) {
-        if (!vectorSampleMean(inF32, errorsF32, maskU8, maskVal, stats)) {
-            psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSampleMean() returned an error.\n");
-            stats->sampleMean = NAN;
-        } else {
-            vectorSampleStdev(inF32, errorsF32, maskU8, maskVal, stats);
-        }
-    }
-
-    // ************************************************************************
-    if (stats->options & (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_ROBUST_QUARTILE |
-                          PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV)) {
-        if (!vectorRobustStats(inF32, errorsF32, maskU8, maskVal, stats)) {
-            psError(PS_ERR_UNKNOWN, false, _("Failed to calculate the specified statistic."));
-            psFree(inF32);
-            psFree(errorsF32);
-            psTrace("psLib.math", 3,"---- %s(NULL) end  ----\n", __func__);
-            return(NULL);
-        }
-    }
-
-    // ************************************************************************
-    if ((stats->options & PS_STAT_CLIPPED_MEAN) || (stats->options & PS_STAT_CLIPPED_STDEV)) {
-        if (!vectorClippedStats(inF32, errorsF32, maskU8, maskVal, stats)) {
-            psError(PS_ERR_UNKNOWN, false, "Failed to calculate clipped statistics for input psVector.\n");
-            stats->clippedMean = NAN;
-            stats->clippedStdev = NAN;
+        if (!vectorSampleStdev(inF32, errorsF32, maskU8, maskVal, stats)) {
+            psError(PS_ERR_UNKNOWN, false, "Failed to calculate sample stdev");
+            status &= false;
         }
     }
@@ -2089 +1545 @@
         if (vectorMinMax(inF32, maskU8, maskVal, stats) == 0) {
             psError(PS_ERR_UNKNOWN, false, "Failed to calculate vector minimum and maximum");
+            status &= false;
+        }
+    }
+
+    // ************************************************************************
+    if (stats->options & (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_ROBUST_QUARTILE)) {
+        if (!vectorRobustStats(inF32, errorsF32, maskU8, maskVal, stats)) {
+            psError(PS_ERR_UNKNOWN, false, _("Failed to calculate robust statistics"));
+            status &= false;
+        }
+    }
+
+    // ************************************************************************
+    if (stats->options & (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV)) {
+        if (!vectorFittedStats(inF32, errorsF32, maskU8, maskVal, stats)) {
+            psError(PS_ERR_UNKNOWN, false, _("Failed to calculate fitted statistics"));
+            status &= false;
+        }
+    }
+
+    // ************************************************************************
+    if (stats->options & (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2)) {
+        if (stats->options & (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV)) {
+            psAbort ("stats", "you may not specify both FITTED_MEAN and FITTED_MEAN_V2");
+        }
+        if (!vectorFittedStats_v2(inF32, errorsF32, maskU8, maskVal, stats)) {
+            psError(PS_ERR_UNKNOWN, false, _("Failed to calculate fitted statistics"));
+            status &= false;
+        }
+    }
+
+    // ************************************************************************
+    if ((stats->options & PS_STAT_CLIPPED_MEAN) || (stats->options & PS_STAT_CLIPPED_STDEV)) {
+        if (!vectorClippedStats(inF32, errorsF32, maskU8, maskVal, stats)) {
+            psError(PS_ERR_UNKNOWN, false, "Failed to calculate clipped statistics\n");
+            status &= false;
         }
     }
@@ -2096 +1588 @@
     psFree(maskU8);
     psTrace("psLib.math", 3,"---- %s() end  ----\n", __func__);
-    return (stats);
+    return (status);
 }
 
@@ -2106 +1598 @@
     }
 
-    READ_STAT("MEAN",     PS_STAT_SAMPLE_MEAN);
-    READ_STAT("STDEV",    PS_STAT_SAMPLE_STDEV);
-    READ_STAT("MEDIAN",   PS_STAT_SAMPLE_MEDIAN);
-    READ_STAT("QUARTILE", PS_STAT_SAMPLE_QUARTILE);
+    READ_STAT("MEAN",       PS_STAT_SAMPLE_MEAN);
+    READ_STAT("STDEV",      PS_STAT_SAMPLE_STDEV);
+    READ_STAT("MEDIAN",     PS_STAT_SAMPLE_MEDIAN);
+    READ_STAT("QUARTILE",   PS_STAT_SAMPLE_QUARTILE);
     READ_STAT("SAMPLE_MEAN",     PS_STAT_SAMPLE_MEAN);
     READ_STAT("SAMPLE_STDEV",    PS_STAT_SAMPLE_STDEV);
@@ -2118 +1610 @@
     READ_STAT("ROBUST_STDEV",    PS_STAT_ROBUST_STDEV);
     READ_STAT("ROBUST_QUARTILE", PS_STAT_ROBUST_QUARTILE);
-    READ_STAT("FITTED",       PS_STAT_FITTED_MEAN);
-    READ_STAT("FITTED_MEAN",  PS_STAT_FITTED_MEAN);
-    READ_STAT("FITTED_STDEV", PS_STAT_ROBUST_STDEV);
-    READ_STAT("CLIPPED",       PS_STAT_CLIPPED_MEAN);
-    READ_STAT("CLIPPED_MEAN",  PS_STAT_CLIPPED_MEAN);
-    READ_STAT("CLIPPED_STDEV", PS_STAT_CLIPPED_STDEV);
+    READ_STAT("FITTED",         PS_STAT_FITTED_MEAN);
+    READ_STAT("FITTED_MEAN",    PS_STAT_FITTED_MEAN);
+    READ_STAT("FITTED_STDEV",   PS_STAT_FITTED_STDEV);
+    READ_STAT("FITTED_V2",       PS_STAT_FITTED_MEAN_V2);
+    READ_STAT("FITTED_MEAN_V2",  PS_STAT_FITTED_MEAN_V2);
+    READ_STAT("FITTED_STDEV_V2", PS_STAT_FITTED_STDEV_V2);
+    READ_STAT("CLIPPED",         PS_STAT_CLIPPED_MEAN);
+    READ_STAT("CLIPPED_MEAN",    PS_STAT_CLIPPED_MEAN);
+    READ_STAT("CLIPPED_STDEV",   PS_STAT_CLIPPED_STDEV);
 
     psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to parse statistic: %s\n", string);
@@ -2146 +1641 @@
     WRITE_STAT("ROBUST_STDEV",    PS_STAT_ROBUST_STDEV);
     WRITE_STAT("ROBUST_QUARTILE", PS_STAT_ROBUST_QUARTILE);
-    WRITE_STAT("FITTED_MEAN",  PS_STAT_FITTED_MEAN);
-    WRITE_STAT("FITTED_STDEV", PS_STAT_ROBUST_STDEV);
-    WRITE_STAT("CLIPPED_MEAN",  PS_STAT_CLIPPED_MEAN);
-    WRITE_STAT("CLIPPED_STDEV", PS_STAT_CLIPPED_STDEV);
+    WRITE_STAT("FITTED_MEAN",     PS_STAT_FITTED_MEAN);
+    WRITE_STAT("FITTED_STDEV",    PS_STAT_FITTED_STDEV);
+    WRITE_STAT("FITTED_MEAN_V2",  PS_STAT_FITTED_MEAN_V2);
+    WRITE_STAT("FITTED_STDEV_V2", PS_STAT_FITTED_STDEV_V2);
+    WRITE_STAT("CLIPPED_MEAN",    PS_STAT_CLIPPED_MEAN);
+    WRITE_STAT("CLIPPED_STDEV",   PS_STAT_CLIPPED_STDEV);
 
     return string;
@@ -2197 +1694 @@
     case PS_STAT_FITTED_MEAN:
     case PS_STAT_FITTED_STDEV:
+    case PS_STAT_FITTED_MEAN_V2:
+    case PS_STAT_FITTED_STDEV_V2:
     case PS_STAT_CLIPPED_MEAN:
     case PS_STAT_CLIPPED_STDEV:
@@ -2227 +1726 @@
     case PS_STAT_FITTED_STDEV:
         return stats->fittedStdev;
+    case PS_STAT_FITTED_MEAN_V2:
+        return stats->fittedMean;
+    case PS_STAT_FITTED_STDEV_V2:
+        return stats->fittedStdev;
     case PS_STAT_CLIPPED_MEAN:
         return stats->clippedMean;
@@ -2246 +1749 @@
     return NAN;
 }
+
+// other private functions used above
+
+static psF32 QuadraticInverse(psF32 a,
+                              psF32 b,
+                              psF32 c,
+                              psF32 y,
+                              psF32 xLo,
+                              psF32 xHi
+                             )
+{
+    psF64 tmp = sqrt((y - c)/a + (b*b)/(4.0 * a * a));
+
+    psF64 x1 = -b/(2.0*a) + tmp;
+    psF64 x2 = -b/(2.0*a) - tmp;
+
+    if (xLo <= x1 && x1 <= xHi) {
+        return x1;
+    }
+    if (xLo <= x2 && x2 <= xHi) {
+        return x2;
+    }
+    return 0.5 * (xLo + xHi);
+}
+
+
+/******************************************************************************
+fitQuadraticSearchForYThenReturnX(*xVec, *yVec, binNum, yVal): A general
+routine which fits a quadratic to three points and returns the x-value
+corresponding to the input y-value.  This routine takes psVectors of x/y pairs
+as input, and fits a quadratic to the 3 points surrounding element binNum in
+the vectors (the midpoint between element i and i+1 is used for x[i]).  It
+then determines for what value x does that quadratic f(x) = yVal (the input
+parameter).
+ 
+*****************************************************************************/
+static psF32 fitQuadraticSearchForYThenReturnX(const psVector *xVec,
+        psVector *yVec,
+        psS32 binNum,
+        psF32 yVal
+                                              )
+{
+    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
+    psTrace("psLib.math", 5, "binNum, yVal is (%d, %f)\n", binNum, yVal);
+    if (psTraceGetLevel("psLib.math") >= 8) {
+        PS_VECTOR_PRINT_F32(xVec);
+        PS_VECTOR_PRINT_F32(yVec);
+    }
+
+    PS_ASSERT_VECTOR_NON_NULL(xVec, NAN);
+    PS_ASSERT_VECTOR_NON_NULL(yVec, NAN);
+    PS_ASSERT_VECTOR_TYPE(xVec, PS_TYPE_F32, NAN);
+    PS_ASSERT_VECTOR_TYPE(yVec, PS_TYPE_F32, NAN);
+    //    PS_ASSERT_VECTORS_SIZE_EQUAL(xVec, yVec, NAN);
+    PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, (int)(xVec->n - 1), NAN);
+    PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, (int)(yVec->n - 1), NAN);
+
+    psVector *x = psVectorAlloc(3, PS_TYPE_F64);
+    psVector *y = psVectorAlloc(3, PS_TYPE_F64);
+    psF32 tmpFloat = 0.0f;
+
+    if ((binNum >= 1) && (binNum < (yVec->n - 2)) && (binNum < (xVec->n - 2))) {
+        // The general case.  We have all three points.
+        x->data.F64[0] = (psF64) (0.5 * (xVec->data.F32[binNum - 1] + xVec->data.F32[binNum]));
+        x->data.F64[1] = (psF64) (0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum+1]));
+        x->data.F64[2] = (psF64) (0.5 * (xVec->data.F32[binNum+1] + xVec->data.F32[binNum+2]));
+        y->data.F64[0] = yVec->data.F32[binNum - 1];
+        y->data.F64[1] = yVec->data.F32[binNum];
+        y->data.F64[2] = yVec->data.F32[binNum + 1];
+        psTrace("psLib.math", 6, "x vec (orig) is (%f %f %f %f)\n", xVec->data.F32[binNum - 1],
+                xVec->data.F32[binNum], xVec->data.F32[binNum+1], xVec->data.F32[binNum+2]);
+        psTrace("psLib.math", 6, "x data is (%f %f %f)\n", x->data.F64[0], x->data.F64[1], x->data.F64[2]);
+        psTrace("psLib.math", 6, "y data is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]);
+
+        //
+        // Ensure that the y value lies within range of the y values.
+        //
+        if (! (((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[2])) ||
+                ((y->data.F64[2] <= yVal) && (yVal <= y->data.F64[0]))) ) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                    _("Specified yVal, %g, is not within y-range, %g to %g."),
+                    (psF64)yVal, y->data.F64[0], y->data.F64[2]);
+        }
+
+        //
+        // Ensure that the y values are monotonic.
+        //
+        if (((y->data.F64[0] < y->data.F64[1]) && !(y->data.F64[1] <= y->data.F64[2])) ||
+                ((y->data.F64[0] > y->data.F64[1]) && !(y->data.F64[1] >= y->data.F64[2]))) {
+            psError(PS_ERR_UNKNOWN, true,
+                    "This routine must be called with monotonically increasing or decreasing data points.\n");
+            psFree(x);
+            psFree(y);
+            psTrace("psLib.math", 5, "---- %s() end ----\n", __func__);
+            return NAN;
+        }
+
+        // Determine the coefficients of the polynomial.
+        psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
+        myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x);
+        if (myPoly == NULL) {
+            psError(PS_ERR_UNEXPECTED_NULL, false,
+                    _("Failed to fit a 1-dimensional polynomial to the three specified data points.  Returning NAN."));
+            psFree(x);
+            psFree(y);
+            psTrace("psLib.math", 5, "---- %s(NAN) end ----\n", __func__);
+            return NAN;
+        }
+        psTrace("psLib.math", 6, "myPoly->coeff[0] is %f\n", myPoly->coeff[0]);
+        psTrace("psLib.math", 6, "myPoly->coeff[1] is %f\n", myPoly->coeff[1]);
+        psTrace("psLib.math", 6, "myPoly->coeff[2] is %f\n", myPoly->coeff[2]);
+        psTrace("psLib.math", 6, "Fitted y vec is (%f %f %f)\n",
+                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[0]),
+                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[1]),
+                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[2]));
+
+        psTrace("psLib.math", 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal);
+        tmpFloat = QuadraticInverse(myPoly->coeff[2], myPoly->coeff[1], myPoly->coeff[0], yVal,
+                                    x->data.F64[0], x->data.F64[2]);
+        psFree(myPoly);
+
+        if (isnan(tmpFloat)) {
+            psError(PS_ERR_UNEXPECTED_NULL,
+                    false, _("Failed to determine the median of the fitted polynomial.  Returning NAN."));
+            psFree(x);
+            psFree(y);
+            psTrace("psLib.math", 5, "---- %s(NAN) end ----\n", __func__);
+            return(NAN);
+        }
+    } else {
+        // These are special cases where the bin is at the beginning or end of the vector.
+        if (binNum == 0) {
+            // We have two points only at the beginning of the vectors x and y.
+            tmpFloat = 0.5 * (xVec->data.F32[binNum] +
+                              xVec->data.F32[binNum + 1]);
+        } else if (binNum == (xVec->n - 1)) {
+            // The special case where we have two points only at the end of
+            // the vectors x and y.
+            // XXX: Is this right?
+            tmpFloat = xVec->data.F32[binNum];
+        } else if (binNum == (xVec->n - 2)) {
+            // XXX: Is this right?
+            tmpFloat = 0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum + 1]);
+        }
+    }
+
+    psTrace("psLib.math", 6, "FIT: return %f\n", tmpFloat);
+    psFree(x);
+    psFree(y);
+
+    psTrace("psLib.math", 5, "---- %s(%f) end ----\n", __func__, tmpFloat);
+    return tmpFloat;
+}
+
+/******************************************************************************
+NOTE: We assume unnormalized gaussians.
+*****************************************************************************/
+static psF32 minimizeLMChi2Gauss1D(psVector *deriv,
+                                   const psVector *params,
+                                   const psVector *coords
+                                  )
+{
+    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
+    PS_ASSERT_VECTOR_NON_NULL(params, NAN);
+    PS_ASSERT_VECTOR_SIZE(params, (long)2, NAN);
+    PS_ASSERT_VECTOR_TYPE(params, PS_TYPE_F32, NAN);
+    PS_ASSERT_VECTOR_NON_NULL(coords, NAN);
+    PS_ASSERT_VECTOR_SIZE(coords, (long)1, NAN);
+    PS_ASSERT_VECTOR_TYPE(coords, PS_TYPE_F32, NAN);
+
+    psF32 x = coords->data.F32[0];
+    psF32 mean = params->data.F32[0];
+    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 = dx * gauss;
+        deriv->data.F32[0] = tmp / var;
+        deriv->data.F32[1] = tmp * dx / (var * var);
+    }
+
+
+    psTrace("psLib.math", 4, "---- %s() end ----\n", __func__);
+    return gauss;
+}

trunk/psLib/src/math/psStats.h

@@ -7 +7 @@
  *  on those data structures.
  *
- *  XXX: The following stats members are never used, or set in this code.
- *      stats->robustN50
- *      stats->clippedNvalues
- *      stats->binsize
- *
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.56 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-08-22 15:02:08 $
+ *  @version $Revision: 1.57 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-12-08 11:38:54 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
  */
+
 #ifndef PS_STATS_H
 #define PS_STATS_H
 
-#include "psVector.h"
+// #include "psVector.h"
 
 /// @addtogroup Stats
@@ -35 +31 @@
  *  @see psStats, psVectorStats, psImageStats
  */
-// XXX: Is PS_STAT_ROBUST_FOR_SAMPLE obsolete?
 typedef enum {
+    PS_STAT_NONE            = 0x000000, ///< Empty set
     PS_STAT_SAMPLE_MEAN     = 0x000001, ///< Sample Mean
     PS_STAT_SAMPLE_MEDIAN   = 0x000002, ///< Sample Median
@@ -51 +47 @@
     PS_STAT_MIN             = 0x001000, ///< Minumum
     PS_STAT_USE_RANGE       = 0x002000, ///< Range
-    PS_STAT_USE_BINSIZE     = 0x004000  ///< Binsize
+    PS_STAT_USE_BINSIZE     = 0x004000, ///< Binsize
+    PS_STAT_FITTED_MEAN_V2  = 0x008000, ///< Fitted Mean
+    PS_STAT_FITTED_STDEV_V2 = 0x010000, ///< Fitted Standard Deviation
 } psStatsOptions;
 
@@ -81 +79 @@
     double binsize;                    ///< binsize for robust fit (input/ouput)
     long nSubsample;                   ///< maxinum number of measurements (input)
-    psStatsOptions options;            ///< bitmask of calculated values
+    psStatsOptions options;            ///< bitmask of values requested
+    psStatsOptions results;            ///< bitmask of values calculated
 }
 psStats;
@@ -89 +88 @@
  *  @return psStats*    the statistical results as specified by stats->options
  */
-psStats* psVectorStats(
+bool psVectorStats(
     psStats* stats,                    ///< stats structure defines stats to be calculated and how
     const psVector* in,                ///< Vector to be analysed.
@@ -117 +116 @@
 );
 
-
-/******************************************************************************
-    Histogram functions and data structures.
- *****************************************************************************/
-
-/** The basic histogram structure which contains bounds and bins.
- *
- *  In this structure, the vector bounds specifies the boundaries of the
- *  histogram bins, and must of type psF32, while nums specifies the number
- *  of entries in the bin, and must of type psU32. The value of bounds.n must
- *  therefore be 1 greater than than nums.n. The two values minNum and maxNum
- *  are the number of data values which fell below the lower limit bound or
- *  above the upper limit bound, respectively.
- */
-typedef struct
-{
-    const psVector* bounds;            ///< Bounds for the bins (type F32)
-    psVector* nums;                    ///< Number in each of the bins (INT)
-    int minNum;                        ///< Number below the minimum
-    int maxNum;                        ///< Number above the maximum
-    bool uniform;                      ///< Is it a uniform distribution?
-}
-psHistogram;
-
-/** Allocator for psHistogram where the bounds of the bins are implicitly
- *  specified through simply specifying an upper and lower limit along with
- *  the size of the bins.
- *
- *  @return psHistogram*    Newly allocated psHistogram
- */
-psHistogram* psHistogramAlloc(
-    float lower,                       ///< Lower limit for the bins
-    float upper,                       ///< Upper limit for the bins
-    int n                              ///< Number of bins
-);
-
-
-/** Checks the type of a particular pointer.
- *
- *  Uses the appropriate deallocation function in psMemBlock to check the ptr datatype.
- *
- *  @return bool:       True if the pointer matches a psHistogram structure, false otherwise.
- */
-bool psMemCheckHistogram(
-    psPtr ptr                          ///< the pointer whose type to check
-);
-
-
-/** Allocator for psHistogram where the bounds of the bins are explicitly
- *  specified.
- *
- *  @return psHistogram*    Newly allocated psHistogram
- */
-psHistogram* psHistogramAllocGeneric(
-    const psVector* bounds             ///< Bounds for the bins
-);
-
-/** Calculate a histogram
- *
- *  The following function populates the histogram bins from the specified
- *  vector (in). It alters and returns the histogram out structure. The input
- *  vector may be of types psU8, psU16, psF32, psF64.
- *
- *  @return psHistogram*   histogram result
- */
-psHistogram* psVectorHistogram(
-    psHistogram* out,                  ///< Histogram data
-    const psVector* values,            ///< Vector to analyse
-    const psVector* errors,            ///< Errors
-    const psVector* mask,              ///< Mask dat for input vector
-    psMaskType maskVal                 ///< Mask value
-);
-
-
 // Get the statistics option from a string
 psStatsOptions psStatsOptionFromString(const char *string);
@@ -207 +132 @@
 
 #endif // #ifndef PS_STATS_H
+
+/*
+ * private stats functions used in psStats.c:
+ *
+ * vectorSampleMean
+   (none)
+ * vectorMinMax
+   (none)
+ * vectorSampleMedian (also yields SAMPLE_QUARTILE)
+   (none) 
+ * vectorSampleStdev
+   (vectorSampleMean)
+ * vectorClippedStats
+   (vectorSampleMedian)
+   (vectorSampleMean (*also subset))
+   (vectorSampleStdev (*also subset))
+ * vectorRobustStats
+   (vectorMinMax (*only subset))
+ * vectorFittedStats
+   (vectorRobustStats)
+ 
+ * private stats functions called by other private stats functions are automatically called by
+ * those functions.  since they set the stats->results flags, they are not called multiple
+ * times.
+ 
+ * the private stats functions do not test for their corresponding stats flags: it is not
+ * necessary to request them if they are called within this function.
+ 
+*/

trunk/psLib/test/math/tap_psStatsTiming.c

@@ -9 +9 @@
 // ok(condition, "condition succeeded");
 // skip_start(condition, Nskip, "Skipping tests because of failure");
+
+# define DTIME(A,B) ((A.tv_sec - B.tv_sec) + 1e-6*(A.tv_usec - B.tv_usec))
+struct timeval start, mark;
 
 int main (void)
@@ -22 +25 @@
         rnd->data.F32[i] = psRandomGaussian (seed);
     }
-    psTimerStart ("test");
-
-    /********** fitted stats ***********/
-    // test FITTED timing
-    {
-        // psMemId id = psMemGetId();
-
-        diag ("check fitted mean speed : 1000 pts / 1000 loops");
+
+    diag ("timing for sample mean");
+    /********** SAMPLE MEAN ***********/
+    // test stat sample mean (no mask, no range)
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN);
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 10000; i++)
+        {
+            psVectorStats (stats, rnd, NULL, NULL, 0);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.1, "sample mean %f (mask: 0, range: 0): %.3f sec", stats->sampleMean, delta);
+        psFree (stats);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+    // test stat sample mean (mask, no range)
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN);
+        psVector *mask = psVectorAlloc (1000, PS_TYPE_U8);
+        psVectorInit (mask, 0);
+        mask->data.U8[100] = 1;
+        mask->data.U8[200] = 1;
+        mask->data.U8[300] = 1;
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 10000; i++)
+        {
+            psVectorStats (stats, rnd, NULL, mask, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.12, "sample mean %f (mask: 1, range: 0): %.3f sec", stats->sampleMean, delta);
+        psFree (stats);
+        psFree (mask);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+    // test stat sample mean (no mask, range)
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_USE_RANGE);
+        stats->min = -10;
+        stats->max = +10;
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 10000; i++)
+        {
+            psVectorStats (stats, rnd, NULL, NULL, 0);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.18, "sample mean %f (mask: 0, range: 1): %.3f sec", stats->sampleMean, delta);
+        psFree (stats);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+    // test stat sample mean (mask, range)
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_USE_RANGE);
+        stats->min = -10;
+        stats->max = +10;
+        psVector *mask = psVectorAlloc (1000, PS_TYPE_U8);
+        psVectorInit (mask, 0);
+        mask->data.U8[100] = 1;
+        mask->data.U8[200] = 1;
+        mask->data.U8[300] = 1;
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 10000; i++)
+        {
+            psVectorStats (stats, rnd, NULL, mask, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.2, "sample mean %f (mask: 1, range: 1): %.3f sec", stats->sampleMean, delta);
+        psFree (stats);
+        psFree (mask);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+    diag ("timing for sample median");
+    /********** SAMPLE MEDIAN ***********/
+    // test stat sample median (no mask, no range)
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 10000; i++)
+        {
+            psVectorStats (stats, rnd, NULL, NULL, 0);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 2.8, "sample median %f (mask: 0, range: 0): %.3f sec", stats->sampleMedian, delta);
+        psFree (stats);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+    // test stat sample median (mask, no range)
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+        psVector *mask = psVectorAlloc (1000, PS_TYPE_U8);
+        psVectorInit (mask, 0);
+        mask->data.U8[100] = 1;
+        mask->data.U8[200] = 1;
+        mask->data.U8[300] = 1;
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 10000; i++)
+        {
+            psVectorStats (stats, rnd, NULL, mask, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 2.8, "sample median %f (mask: 1, range: 0): %.3f sec", stats->sampleMedian, delta);
+        psFree (stats);
+        psFree (mask);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+    // test stat sample median (no mask, range)
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_USE_RANGE);
+        stats->min = -10;
+        stats->max = +10;
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 10000; i++)
+        {
+            psVectorStats (stats, rnd, NULL, NULL, 0);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 2.8, "sample median %f (mask: 0, range: 1): %.3f sec", stats->sampleMedian, delta);
+        psFree (stats);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+    // test stat sample median (mask, range)
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_USE_RANGE);
+        stats->min = -10;
+        stats->max = +10;
+        psVector *mask = psVectorAlloc (1000, PS_TYPE_U8);
+        psVectorInit (mask, 0);
+        mask->data.U8[100] = 1;
+        mask->data.U8[200] = 1;
+        mask->data.U8[300] = 1;
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 10000; i++)
+        {
+            psVectorStats (stats, rnd, NULL, mask, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 2.8, "sample median %f (mask: 1, range: 1): %.3f sec", stats->sampleMedian, delta);
+        psFree (stats);
+        psFree (mask);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+    diag ("timing for sample stdev");
+    /********** SAMPLE STDEV ***********/
+    // test stat sample stdev (no mask, no range)
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_STDEV);
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 10000; i++)
+        {
+            psVectorStats (stats, rnd, NULL, NULL, 0);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.2, "sample stdev %f (mask: 0, range: 0): %.3f sec", stats->sampleStdev, delta);
+        psFree (stats);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+    // test stat sample stdev (mask, no range)
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_STDEV);
+        psVector *mask = psVectorAlloc (1000, PS_TYPE_U8);
+        psVectorInit (mask, 0);
+        mask->data.U8[100] = 1;
+        mask->data.U8[200] = 1;
+        mask->data.U8[300] = 1;
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 10000; i++)
+        {
+            psVectorStats (stats, rnd, NULL, mask, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.27, "sample stdev %f (mask: 1, range: 0): %.3f sec", stats->sampleStdev, delta);
+        psFree (stats);
+        psFree (mask);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+    // test stat sample stdev (no mask, range)
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_STDEV | PS_STAT_USE_RANGE);
+        stats->min = -10;
+        stats->max = +10;
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 10000; i++)
+        {
+            psVectorStats (stats, rnd, NULL, NULL, 0);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.36, "sample stdev %f (mask: 0, range: 1): %.3f sec", stats->sampleStdev, delta);
+        psFree (stats);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+    // test stat sample stdev (mask, range)
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_STDEV | PS_STAT_USE_RANGE);
+        stats->min = -10;
+        stats->max = +10;
+        psVector *mask = psVectorAlloc (1000, PS_TYPE_U8);
+        psVectorInit (mask, 0);
+        mask->data.U8[100] = 1;
+        mask->data.U8[200] = 1;
+        mask->data.U8[300] = 1;
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 10000; i++)
+        {
+            psVectorStats (stats, rnd, NULL, mask, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.42, "sample stdev %f (mask: 1, range: 1): %.3f sec", stats->sampleStdev, delta);
+        psFree (stats);
+        psFree (mask);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+    diag ("timing for sample min,max");
+    /*************** MIN,MAX ******************/
+    // test stat min,max (no mask, no range)
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 10000; i++)
+        {
+            psVectorStats (stats, rnd, NULL, NULL, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.17, "sample min,max %f,%f (mask: 0, range: 0): %.3f sec", stats->min, stats->max, delta);
+        psFree (stats);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+    // test stat min,max (no mask, no range)
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
+        psVector *mask = psVectorAlloc (1000, PS_TYPE_U8);
+        psVectorInit (mask, 0);
+        mask->data.U8[100] = 1;
+        mask->data.U8[200] = 1;
+        mask->data.U8[300] = 1;
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 10000; i++)
+        {
+            psVectorStats (stats, rnd, NULL, mask, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.18, "sample min,max %f,%f (mask: 1, range: 0): %.3f sec", stats->min, stats->max, delta);
+        psFree (stats);
+        psFree (mask);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+    // test stat min,max (no mask, no range)
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX | PS_STAT_USE_RANGE);
+        stats->min = -10;
+        stats->max = +10;
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 10000; i++)
+        {
+            psVectorStats (stats, rnd, NULL, NULL, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.22, "sample min,max %f,%f (mask: 0, range: 1): %.3f sec", stats->min, stats->max, delta);
+        psFree (stats);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+    // test stat min,max (no mask, no range)
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX | PS_STAT_USE_RANGE);
+        stats->min = -10;
+        stats->max = +10;
+        psVector *mask = psVectorAlloc (1000, PS_TYPE_U8);
+        psVectorInit (mask, 0);
+        mask->data.U8[100] = 1;
+        mask->data.U8[200] = 1;
+        mask->data.U8[300] = 1;
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 10000; i++)
+        {
+            psVectorStats (stats, rnd, NULL, mask, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.26, "sample min,max %f,%f (mask: 1, range: 1): %.3f sec", stats->min, stats->max, delta);
+        psFree (stats);
+        psFree (mask);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+    diag ("timing for clipped stats");
+    /********** CLIPPED STATS ***********/
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
+        psVector *rnd2 = psVectorAlloc (1000, PS_TYPE_F32);
+        for (int i = 0; i < rnd2->n; i++)
+        {
+            rnd2->data.F32[i] = psRandomGaussian (seed);
+        }
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 1000; i++)
+        {
+            psVectorStats (stats, rnd2, NULL, NULL, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.3, "clipped mean %f, stdev %f (mask: 0, range: 0): %.3f sec (1000 pts / 1000 loops)", stats->clippedMean, stats->clippedStdev, delta);
+        psFree (stats);
+        psFree (rnd2);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
+        psVector *rnd2 = psVectorAlloc (3000, PS_TYPE_F32);
+        for (int i = 0; i < rnd2->n; i++)
+        {
+            rnd2->data.F32[i] = psRandomGaussian (seed);
+        }
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 1000; i++)
+        {
+            psVectorStats (stats, rnd2, NULL, NULL, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.5, "clipped mean %f, stdev %f (mask: 0, range: 0): %.3f sec (3000 pts / 1000 loops)", stats->clippedMean, stats->clippedStdev, delta);
+        psFree (stats);
+        psFree (rnd2);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
+        psVector *rnd2 = psVectorAlloc (10000, PS_TYPE_F32);
+        for (int i = 0; i < rnd2->n; i++)
+        {
+            rnd2->data.F32[i] = psRandomGaussian (seed);
+        }
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 1000; i++)
+        {
+            psVectorStats (stats, rnd2, NULL, NULL, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 1.2, "clipped mean %f, stdev %f (mask: 0, range: 0): %.3f sec (10000 pts / 1000 loops)", stats->clippedMean, stats->clippedStdev, delta);
+        psFree (stats);
+        psFree (rnd2);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+    diag ("timing for robust stats");
+    /********** ROBUST STATS ***********/
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_ROBUST_QUARTILE);
+        psVector *rnd2 = psVectorAlloc (1000, PS_TYPE_F32);
+        for (int i = 0; i < rnd2->n; i++)
+        {
+            rnd2->data.F32[i] = psRandomGaussian (seed);
+        }
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 1000; i++)
+        {
+            psVectorStats (stats, rnd2, NULL, NULL, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.3, "robust mean %f, stdev %f (mask: 0, range: 0): %.3f sec (1000 pts / 1000 loops)", stats->robustMedian, stats->robustStdev, delta);
+        psFree (stats);
+        psFree (rnd2);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_ROBUST_QUARTILE);
+        psVector *rnd2 = psVectorAlloc (3000, PS_TYPE_F32);
+        for (int i = 0; i < rnd2->n; i++)
+        {
+            rnd2->data.F32[i] = psRandomGaussian (seed);
+        }
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 1000; i++)
+        {
+            psVectorStats (stats, rnd2, NULL, NULL, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.5, "robust mean %f, stdev %f (mask: 0, range: 0): %.3f sec (3000 pts / 1000 loops)", stats->robustMedian, stats->robustStdev, delta);
+        psFree (stats);
+        psFree (rnd2);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_ROBUST_QUARTILE);
+        psVector *rnd2 = psVectorAlloc (10000, PS_TYPE_F32);
+        for (int i = 0; i < rnd2->n; i++)
+        {
+            rnd2->data.F32[i] = psRandomGaussian (seed);
+        }
+
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 1000; i++)
+        {
+            psVectorStats (stats, rnd2, NULL, NULL, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 1.2, "robust mean %f, stdev %f (mask: 0, range: 0): %.3f sec (10000 pts / 1000 loops)", stats->robustMedian, stats->robustStdev, delta);
+        psFree (stats);
+        psFree (rnd2);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+    diag ("timing for fitted stats");
+    /********** FITTED TIMING ***********/
+    {
+        psMemId id = psMemGetId();
 
         psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV);
@@ -38 +554 @@
         }
 
-        psTimerClear ("test");
-        for (int i = 0; i < 1000; i++)
-        {
-            psVectorStats (stats, rnd2, NULL, NULL, 1);
-        }
-        psF64 delta = psTimerMark("test");
-        ok (delta < 0.5, "fitted mean %f, stdev %f (mask: 0, range: 0): %.3f sec", stats->fittedMean, stats->fittedStdev, delta);
-        psFree (stats);
-
-        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-    // test FITTED timing
-    {
-        // psMemId id = psMemGetId();
-
-        diag ("check fitted mean speed : 3000 pts / 1000 loops");
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 1000; i++)
+        {
+            psVectorStats (stats, rnd2, NULL, NULL, 1);
+
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.7, "fitted mean %f, stdev %f (mask: 0, range: 0): %.3f sec (1000 pts / 1000 loops)", stats->fittedMean, stats->fittedStdev, delta);
+        psFree (stats);
+        psFree (rnd2);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+    {
+        psMemId id = psMemGetId();
 
         psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV);
@@ -62 +578 @@
         }
 
-        psTimerClear ("test");
-        for (int i = 0; i < 1000; i++)
-        {
-            psVectorStats (stats, rnd2, NULL, NULL, 1);
-        }
-        psF64 delta = psTimerMark("test");
-        ok (delta < 0.8, "fitted mean %f, stdev %f (mask: 0, range: 0): %.3f sec", stats->fittedMean, stats->fittedStdev, delta);
-        psFree (stats);
-
-        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-    // test FITTED timing
-    {
-        // psMemId id = psMemGetId();
-
-        diag ("check fitted mean speed : 10000 pts / 1000 loops");
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 1000; i++)
+        {
+            psVectorStats (stats, rnd2, NULL, NULL, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.8, "fitted mean %f, stdev %f (mask: 0, range: 0): %.3f sec (3000 pts / 1000 loops)", stats->fittedMean, stats->fittedStdev, delta);
+        psFree (stats);
+        psFree (rnd2);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+    {
+        psMemId id = psMemGetId();
 
         psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV);
@@ -86 +601 @@
         }
 
-        psTimerClear ("test");
-        for (int i = 0; i < 1000; i++)
-        {
-            psVectorStats (stats, rnd2, NULL, NULL, 1);
-        }
-        psF64 delta = psTimerMark("test");
-        ok (delta < 2.2, "fitted mean %f, stdev %f (mask: 0, range: 0): %.3f sec", stats->fittedMean, stats->fittedStdev, delta);
-        psFree (stats);
-
-        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-    // test stat min,max (no mask, no range)
-    {
-        // psMemId id = psMemGetId();
-
-        diag ("check fitted mean speed (and value)");
-
-        psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV | PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
-        psVector *sample = psVectorAlloc (1000, PS_TYPE_F32);
-        psVector *fitted = psVectorAlloc (1000, PS_TYPE_F32);
-
-        for (int i = 0; i < 1000; i++)
-        {
-            // generate a new sample
-            for (int j = 0; j < rnd->n; j++) {
-                rnd->data.F32[j] = psRandomGaussian (seed);
-            }
-            // measure the stats
-            psVectorStats (stats, rnd, NULL, NULL, 1);
-            sample->data.F32[i] = stats->sampleMean;
-            fitted->data.F32[i] = stats->fittedMean;
-        }
-        psFree (stats);
-
-        stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
-        psVectorStats (stats, sample, NULL, NULL, 1);
-        ok (stats->sampleStdev < 2/sqrt(1000), "sample mean %f, stdev %f", stats->sampleMean, stats->sampleStdev);
-        psVectorStats (stats, fitted, NULL, NULL, 1);
-        ok (stats->sampleStdev < 2/sqrt(1000), "fitted mean %f, stdev %f", stats->sampleMean, stats->sampleStdev);
-        psFree (stats);
-
-        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-
-    /********** robust stats ***********/
-    // test ROBUST timing
-    {
-        // psMemId id = psMemGetId();
-
-        diag ("check robust mean speed : 1000 pts / 1000 loops");
-
-        psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_ROBUST_QUARTILE);
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 1000; i++)
+        {
+            psVectorStats (stats, rnd2, NULL, NULL, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 2.2, "fitted mean %f, stdev %f (mask: 0, range: 0): %.3f sec (10000 pts / 1000 loops)", stats->fittedMean, stats->fittedStdev, delta);
+        psFree (stats);
+        psFree (rnd2);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+    diag ("timing for fitted (v2) stats");
+    /********** FITTED (v2) TIMING ***********/
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2);
         psVector *rnd2 = psVectorAlloc (1000, PS_TYPE_F32);
         for (int i = 0; i < rnd2->n; i++)
@@ -144 +627 @@
         }
 
-        psTimerClear ("test");
-        for (int i = 0; i < 1000; i++)
-        {
-            psVectorStats (stats, rnd2, NULL, NULL, 1);
-        }
-        psF64 delta = psTimerMark("test");
-        ok (delta < 0.2, "robust mean %f, stdev %f (mask: 0, range: 0): %.3f sec", stats->robustMedian, stats->robustStdev, delta);
-        psFree (stats);
-
-        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-    // test ROBUST timing
-    {
-        // psMemId id = psMemGetId();
-
-        diag ("check robust mean speed : 3000 pts / 1000 loops");
-
-        psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_ROBUST_QUARTILE);
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 1000; i++)
+        {
+            psVectorStats (stats, rnd2, NULL, NULL, 1);
+
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.7, "fitted mean %f, stdev %f (mask: 0, range: 0): %.3f sec (1000 pts / 1000 loops)", stats->fittedMean, stats->fittedStdev, delta);
+        psFree (stats);
+        psFree (rnd2);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2);
         psVector *rnd2 = psVectorAlloc (3000, PS_TYPE_F32);
         for (int i = 0; i < rnd2->n; i++)
@@ -168 +651 @@
         }
 
-        psTimerClear ("test");
-        for (int i = 0; i < 1000; i++)
-        {
-            psVectorStats (stats, rnd2, NULL, NULL, 1);
-        }
-        psF64 delta = psTimerMark("test");
-        ok (delta < 0.5, "robust mean %f, stdev %f (mask: 0, range: 0): %.3f sec", stats->robustMedian, stats->robustStdev, delta);
-        psFree (stats);
-
-        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-    // test ROBUST timing
-    {
-        // psMemId id = psMemGetId();
-
-        diag ("check robust mean speed : 10000 pts / 1000 loops");
-
-        psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_ROBUST_QUARTILE);
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 1000; i++)
+        {
+            psVectorStats (stats, rnd2, NULL, NULL, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 0.8, "fitted mean %f, stdev %f (mask: 0, range: 0): %.3f sec (3000 pts / 1000 loops)", stats->fittedMean, stats->fittedStdev, delta);
+        psFree (stats);
+        psFree (rnd2);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2);
         psVector *rnd2 = psVectorAlloc (10000, PS_TYPE_F32);
         for (int i = 0; i < rnd2->n; i++)
@@ -192 +674 @@
         }
 
-        psTimerClear ("test");
-        for (int i = 0; i < 1000; i++)
-        {
-            psVectorStats (stats, rnd2, NULL, NULL, 1);
-        }
-        psF64 delta = psTimerMark("test");
-        ok (delta < 1.2, "robust mean %f, stdev %f (mask: 0, range: 0): %.3f sec", stats->robustMedian, stats->robustStdev, delta);
-        psFree (stats);
-
-        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-    // test stat min,max (no mask, no range)
-    {
-        // psMemId id = psMemGetId();
-
-        diag ("check robust mean speed (and value)");
-
-        psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
+        gettimeofday (&start, NULL);
+        for (int i = 0; i < 1000; i++)
+        {
+            psVectorStats (stats, rnd2, NULL, NULL, 1);
+        }
+        gettimeofday (&mark, NULL);
+        psF64 delta = DTIME(mark, start);
+        ok (delta < 2.2, "fitted mean %f, stdev %f (mask: 0, range: 0): %.3f sec (10000 pts / 1000 loops)", stats->fittedMean, stats->fittedStdev, delta);
+        psFree (stats);
+        psFree (rnd2);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+    diag ("compare sample, robust, and fitted mean and stdev to theoretical");
+    // compare SAMPLE, FITTED, ROBUST mean to theoretical
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV);
         psVector *sample = psVectorAlloc (1000, PS_TYPE_F32);
         psVector *robust = psVectorAlloc (1000, PS_TYPE_F32);
+        psVector *fitted = psVectorAlloc (1000, PS_TYPE_F32);
 
         for (int i = 0; i < 1000; i++)
@@ -223 +708 @@
             sample->data.F32[i] = stats->sampleMean;
             robust->data.F32[i] = stats->robustMedian;
+            fitted->data.F32[i] = stats->fittedMean;
         }
         psFree (stats);
@@ -228 +714 @@
         stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
         psVectorStats (stats, sample, NULL, NULL, 1);
-        ok (stats->sampleStdev < 2/sqrt(1000), "sample mean %f, stdev %f", stats->sampleMean, stats->sampleStdev);
+        ok (stats->sampleStdev < 2/sqrt(1000), "sample mean %f, stdev %f (1000 tries)", stats->sampleMean, stats->sampleStdev);
         psVectorStats (stats, robust, NULL, NULL, 1);
-        ok (stats->sampleStdev < 2/sqrt(1000), "robust mean %f, stdev %f", stats->sampleMean, stats->sampleStdev);
-        psFree (stats);
-
-        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-
-    /**********  stats ***********/
-    // test stat sample mean (no mask, no range)
-    {
-        // psMemId id = psMemGetId();
-
-        diag ("check sample mean speed (and value)");
-
-        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN);
-
-        psTimerClear ("test");
-        for (int i = 0; i < 10000; i++)
-        {
-            psVectorStats (stats, rnd, NULL, NULL, 0);
-        }
-        psF64 delta = psTimerMark("test");
-        ok (delta < 0.1, "sample mean %f (mask: 0, range: 0): %.3f sec", stats->sampleMean, delta);
-        psFree (stats);
-
-        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-
-    // test stat sample mean (mask, no range)
-    {
-        // psMemId id = psMemGetId();
-
-        diag ("check sample mean speed (and value)");
-
-        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN);
-        psVector *mask = psVectorAlloc (1000, PS_TYPE_U8);
-        psVectorInit (mask, 0);
-        mask->data.U8[100] = 1;
-        mask->data.U8[200] = 1;
-        mask->data.U8[300] = 1;
-
-        psTimerClear ("test");
-        for (int i = 0; i < 10000; i++)
-        {
-            psVectorStats (stats, rnd, NULL, mask, 1);
-        }
-        psF64 delta = psTimerMark("test");
-        ok (delta < 0.12, "sample mean %f (mask: 1, range: 0): %.3f sec", stats->sampleMean, delta);
-        psFree (stats);
-
-        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-
-    // test stat sample mean (no mask, range)
-    {
-        // psMemId id = psMemGetId();
-
-        diag ("check sample mean speed (and value)");
-
-        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_USE_RANGE);
-        stats->min = -10;
-        stats->max = +10;
-
-        psTimerClear ("test");
-        for (int i = 0; i < 10000; i++)
-        {
-            psVectorStats (stats, rnd, NULL, NULL, 0);
-        }
-        psF64 delta = psTimerMark("test");
-        ok (delta < 0.18, "sample mean %f (mask: 0, range: 1): %.3f sec", stats->sampleMean, delta);
-        psFree (stats);
-
-        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-
-    // test stat sample mean (mask, range)
-    {
-        // psMemId id = psMemGetId();
-
-        diag ("check sample mean speed (and value)");
-
-        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_USE_RANGE);
-        stats->min = -10;
-        stats->max = +10;
-        psVector *mask = psVectorAlloc (1000, PS_TYPE_U8);
-        psVectorInit (mask, 0);
-        mask->data.U8[100] = 1;
-        mask->data.U8[200] = 1;
-        mask->data.U8[300] = 1;
-
-        psTimerClear ("test");
-        for (int i = 0; i < 10000; i++)
-        {
-            psVectorStats (stats, rnd, NULL, mask, 1);
-        }
-        psF64 delta = psTimerMark("test");
-        ok (delta < 0.2, "sample mean %f (mask: 1, range: 1): %.3f sec", stats->sampleMean, delta);
-        psFree (stats);
-
-        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-
-    /*************** min,max ******************/
-    // test stat min,max (no mask, no range)
-    {
-        // psMemId id = psMemGetId();
-
-        diag ("check min,max speed (no mask, no range)");
-
-        psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
-        psTimerClear ("test");
-        for (int i = 0; i < 10000; i++)
-        {
+        ok (stats->sampleStdev < 2/sqrt(1000), "robust mean %f, stdev %f (1000 tries)", stats->sampleMean, stats->sampleStdev);
+        psVectorStats (stats, fitted, NULL, NULL, 1);
+        ok (stats->sampleStdev < 2/sqrt(1000), "fitted mean %f, stdev %f (1000 tries)", stats->sampleMean, stats->sampleStdev);
+        psFree (stats);
+        psFree (sample);
+        psFree (robust);
+        psFree (fitted);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+    diag ("compare sample, robust, and fitted mean and stdev to theoretical");
+    // compare SAMPLE, FITTED_V2, ROBUST mean to theoretical
+    {
+        psMemId id = psMemGetId();
+
+        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2);
+        psVector *sample = psVectorAlloc (1000, PS_TYPE_F32);
+        psVector *robust = psVectorAlloc (1000, PS_TYPE_F32);
+        psVector *fitted = psVectorAlloc (1000, PS_TYPE_F32);
+
+        for (int i = 0; i < 1000; i++)
+        {
+            // generate a new sample
+            for (int j = 0; j < rnd->n; j++) {
+                rnd->data.F32[j] = psRandomGaussian (seed);
+            }
+            // measure the stats
             psVectorStats (stats, rnd, NULL, NULL, 1);
-        }
-        psF64 delta = psTimerMark("test");
-        ok (delta < 0.17, "sample min,max %f,%f (mask: 0, range: 0): %.3f sec", stats->min, stats->max, delta);
-        psFree (stats);
-
-        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-    // test stat min,max (no mask, no range)
-    {
-        // psMemId id = psMemGetId();
-
-        diag ("check min,max speed (mask, no range)");
-
-        psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
-        psVector *mask = psVectorAlloc (1000, PS_TYPE_U8);
-        psVectorInit (mask, 0);
-        mask->data.U8[100] = 1;
-        mask->data.U8[200] = 1;
-        mask->data.U8[300] = 1;
-
-        psTimerClear ("test");
-        for (int i = 0; i < 10000; i++)
-        {
-            psVectorStats (stats, rnd, NULL, mask, 1);
-        }
-        psF64 delta = psTimerMark("test");
-        ok (delta < 0.18, "sample min,max %f,%f (mask: 1, range: 0): %.3f sec", stats->min, stats->max, delta);
-        psFree (stats);
-
-        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-    // test stat min,max (no mask, no range)
-    {
-        // psMemId id = psMemGetId();
-
-        diag ("check min,max speed (no mask, range)");
-
-        psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX | PS_STAT_USE_RANGE);
-        stats->min = -10;
-        stats->max = +10;
-
-        psTimerClear ("test");
-        for (int i = 0; i < 10000; i++)
-        {
-            psVectorStats (stats, rnd, NULL, NULL, 1);
-        }
-        psF64 delta = psTimerMark("test");
-        ok (delta < 0.22, "sample min,max %f,%f (mask: 0, range: 1): %.3f sec", stats->min, stats->max, delta);
-        psFree (stats);
-
-        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
-    }
-    // test stat min,max (no mask, no range)
-    {
-        // psMemId id = psMemGetId();
-
-        diag ("check min,max speed (mask, range)");
-
-        psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX | PS_STAT_USE_RANGE);
-        stats->min = -10;
-        stats->max = +10;
-        psVector *mask = psVectorAlloc (1000, PS_TYPE_U8);
-        psVectorInit (mask, 0);
-        mask->data.U8[100] = 1;
-        mask->data.U8[200] = 1;
-        mask->data.U8[300] = 1;
-
-        psTimerClear ("test");
-        for (int i = 0; i < 10000; i++)
-        {
-            psVectorStats (stats, rnd, NULL, mask, 1);
-        }
-        psF64 delta = psTimerMark("test");
-        ok (delta < 0.26, "sample min,max %f,%f (mask: 1, range: 1): %.3f sec", stats->min, stats->max, delta);
-        psFree (stats);
-
-        // ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+            sample->data.F32[i] = stats->sampleMean;
+            robust->data.F32[i] = stats->robustMedian;
+            fitted->data.F32[i] = stats->fittedMean;
+        }
+        psFree (stats);
+
+        stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
+        psVectorStats (stats, sample, NULL, NULL, 1);
+        ok (stats->sampleStdev < 2/sqrt(1000), "sample mean %f, stdev %f (1000 tries)", stats->sampleMean, stats->sampleStdev);
+        psVectorStats (stats, robust, NULL, NULL, 1);
+        ok (stats->sampleStdev < 2/sqrt(1000), "robust mean %f, stdev %f (1000 tries)", stats->sampleMean, stats->sampleStdev);
+        psVectorStats (stats, fitted, NULL, NULL, 1);
+        ok (stats->sampleStdev < 2/sqrt(1000), "fitted mean %f, stdev %f (1000 tries)", stats->sampleMean, stats->sampleStdev);
+        psFree (stats);
+        psFree (sample);
+        psFree (robust);
+        psFree (fitted);
+
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
     }
 

trunk/psLib/test/math/tap_psStats_Sample_01.c

@@ -215 +215 @@
                          };
 
+static float yraw_03[] = {
+                             183.6000061, 171.3600006, 182.5800018, 178.5, 176.4600067,
+                             176.4600067, 171.3600006, 178.5, 183.6000061, 175.4400024,
+                             159.1199951, 180.5399933, 183.6000061, 187.6799927, 175.4400024,
+                             181.5599976, 174.4199982, 179.5200043, 180.5399933, 187.6799927,
+                             183.6000061, 171.3600006, 177.4799957, 192.7799988, 184.6199951,
+                             176.4600067, 178.5, 182.5800018, 173.3999939, 191.7599945,
+                             202.9799957, 167.2799988, 181.5599976, 185.6399994, 173.3999939,
+                             159.1199951, 178.5, 188.6999969, 179.5200043, 178.5,
+                             184.6199951, 175.4400024, 186.6600037, 179.5200043, 183.6000061,
+                             182.5800018, 173.3999939, 181.5599976, 184.6199951, 177.4799957,
+                             162.1799927, 186.6600037, 163.1999969, 199.9199982, 173.3999939,
+                             170.3399963, 184.6199951, 190.7400055, 173.3999939, 178.5,
+                             189.7200012, 192.7799988, 166.2599945, 187.6799927, 174.4199982,
+                             176.4600067, 181.5599976, 199.9199982, 193.8000031, 170.3399963,
+                             182.5800018, 199.9199982, 183.6000061, 160.1399994, 184.6199951,
+                             195.8399963, 173.3999939, 189.7200012, 192.7799988, 162.1799927,
+                             174.4199982, 190.7400055, 175.4400024, 192.7799988, 181.5599976,
+                             185.6399994, 184.6199951, 166.2599945, 189.7200012, 192.7799988,
+                             183.6000061, 190.7400055, 188.6999969, 174.4199982, 184.6199951,
+                             186.6600037, 182.5800018, 173.3999939, 171.3600006, 166.2599945,
+                             166.2599945, 176.4600067, 147.8999939, 187.6799927, 166.2599945,
+                             173.3999939, 179.5200043, 166.2599945, 169.3200073, 172.3800049,
+                             174.4199982, 195.8399963, 192.7799988, 177.4799957, 181.5599976,
+                             163.1999969, 171.3600006, 194.8200073, 159.1199951, 173.3999939,
+                             182.5800018, 178.5, 179.5200043, 184.6199951, 179.5200043,
+                             164.2200012, 191.7599945, 177.4799957, 186.6600037, 168.3000031,
+                             168.3000031, 168.3000031, 180.5399933, 179.5200043, 175.4400024,
+                             174.4199982, 184.6199951, 186.6600037, 164.2200012, 169.3200073,
+                             187.6799927, 189.7200012, 166.2599945, 169.3200073, 187.6799927,
+                             172.3800049, 204, 194.8200073, 177.4799957, 188.6999969,
+                             171.3600006, 184.6199951, 182.5800018, 184.6199951, 176.4600067,
+                             186.6600037, 170.3399963, 170.3399963, 182.5800018, 184.6199951,
+                             168.3000031, 175.4400024, 176.4600067, 164.2200012, 183.6000061,
+                             161.1600037, 175.4400024, 167.2799988, 179.5200043, 176.4600067,
+                             175.4400024, 179.5200043, 158.1000061, 168.3000031, 168.3000031,
+                             179.5200043, 192.7799988, 190.7400055, 178.5, 173.3999939,
+                             180.5399933, 179.5200043, 191.7599945, 184.6199951, 177.4799957,
+                             182.5800018, 188.6999969, 169.3200073, 162.1799927, 172.3800049,
+                             160.1399994, 177.4799957, 176.4600067, 180.5399933, 170.3399963,
+                             168.3000031, 171.3600006, 180.5399933, 162.1799927, 163.1999969,
+                             170.3399963, 180.5399933, 164.2200012, 183.6000061, 178.5,
+                             193.8000031, 177.4799957, 178.5, 183.6000061, 192.7799988,
+                             157.0800018, 176.4600067, 190.7400055, 179.5200043, 157.0800018,
+                             178.5, 165.2400055, 195.8399963, 177.4799957, 166.2599945,
+                             172.3800049, 182.5800018, 175.4400024, 187.6799927, 164.2200012,
+                             178.5, 173.3999939, 173.3999939, 173.3999939, 172.3800049,
+                             183.6000061, 181.5599976, 163.1999969, 179.5200043, 193.8000031,
+                             172.3800049, 158.1000061, 171.3600006, 165.2400055, 183.6000061,
+                             187.6799927, 158.1000061, 171.3600006, 168.3000031, 166.2599945,
+                             166.2599945, 178.5, 179.5200043, 165.2400055, 176.4600067,
+                             160.1399994, 193.8000031, 170.3399963, 162.1799927, 180.5399933,
+                             159.1199951, 189.7200012, 160.1399994, 196.8600006, 173.3999939,
+                             170.3399963, 165.2400055, 193.8000031, 192.7799988, 184.6199951,
+                             151.9799957, 167.2799988, 162.1799927, 183.6000061, 166.2599945,
+                             177.4799957, 168.3000031, 191.7599945, 181.5599976, 166.2599945,
+                             171.3600006, 188.6999969, 172.3800049, 195.8399963, 189.7200012,
+                             177.4799957, 167.2799988, 191.7599945, 169.3200073, 181.5599976,
+                             180.5399933, 162.1799927, 171.3600006, 162.1799927, 171.3600006,
+                             177.4799957, 167.2799988, 180.5399933, 182.5800018, 187.6799927,
+                             175.4400024, 193.8000031, 174.4199982, 183.6000061, 179.5200043,
+                             164.2200012, 162.1799927, 186.6600037, 189.7200012, 186.6600037,
+                             171.3600006, 187.6799927, 189.7200012, 175.4400024, 184.6199951,
+                             172.3800049, 177.4799957, 183.6000061, 176.4600067, 171.3600006,
+                             172.3800049, 181.5599976, 175.4400024, 170.3399963, 196.8600006,
+                             201.9600067, 180.5399933, 176.4600067, 183.6000061, 175.4400024,
+                             184.6199951, 163.1999969, 177.4799957, 176.4600067, 175.4400024,
+                             192.7799988, 167.2799988, 176.4600067, 184.6199951, 184.6199951,
+                             177.4799957, 185.6399994, 196.8600006, 186.6600037, 184.6199951,
+                             169.3200073, 164.2200012, 185.6399994, 168.3000031, 186.6600037,
+                             180.5399933, 177.4799957, 177.4799957, 186.6600037, 191.7599945,
+                             177.4799957, 181.5599976, 166.2599945, 179.5200043, 184.6199951,
+                             178.5, 171.3600006, 184.6199951, 192.7799988, 182.5800018,
+                             186.6600037, 175.4400024, 182.5800018, 173.3999939, 186.6600037,
+                             162.1799927, 175.4400024, 196.8600006, 166.2599945, 189.7200012,
+                             189.7200012, 182.5800018, 172.3800049, 165.2400055, 191.7599945,
+                             174.4199982, 173.3999939, 178.5, 171.3600006, 157.0800018,
+                             182.5800018, 161.1600037, 195.8399963, 186.6600037, 167.2799988,
+                             166.2599945, 182.5800018, 178.5, 176.4600067, 181.5599976,
+                             184.6199951, 183.6000061, 179.5200043, 174.4199982, 167.2799988,
+                             187.6799927, 176.4600067, 165.2400055, 179.5200043, 157.0800018,
+                             171.3600006, 170.3399963, 175.4400024, 161.1600037, 185.6399994,
+                             169.3200073, 192.7799988, 175.4400024, 172.3800049, 180.5399933,
+                             183.6000061, 174.4199982, 176.4600067, 164.2200012, 183.6000061,
+                             179.5200043, 165.2400055, 169.3200073, 172.3800049, 149.9400024,
+                             175.4400024, 188.6999969, 190.7400055, 171.3600006, 172.3800049,
+                             183.6000061, 178.5, 165.2400055, 176.4600067, 177.4799957,
+                             188.6999969, 192.7799988, 183.6000061, 163.1999969, 186.6600037,
+                             183.6000061, 160.1399994, 167.2799988, 172.3800049, 179.5200043,
+                             189.7200012, 172.3800049, 177.4799957, 161.1600037, 174.4199982,
+                             190.7400055, 181.5599976, 187.6799927, 176.4600067, 176.4600067,
+                             183.6000061, 176.4600067, 188.6999969, 174.4199982, 170.3399963,
+                             185.6399994, 177.4799957, 172.3800049, 179.5200043, 175.4400024,
+                             170.3399963, 169.3200073, 176.4600067, 177.4799957, 169.3200073,
+                             199.9199982, 171.3600006, 194.8200073, 188.6999969, 193.8000031,
+                             182.5800018, 171.3600006, 177.4799957, 175.4400024, 172.3800049,
+                             166.2599945, 183.6000061, 157.0800018, 177.4799957, 193.8000031,
+                             168.3000031, 175.4400024, 175.4400024, 170.3399963, 191.7599945,
+                             189.7200012, 182.5800018, 177.4799957, 157.0800018, 174.4199982,
+                             189.7200012, 162.1799927, 184.6199951, 164.2200012, 157.0800018,
+                             197.8800049, 175.4400024, 184.6199951, 202.9799957, 190.7400055,
+                             171.3600006, 160.1399994, 162.1799927, 176.4600067, 180.5399933,
+                             206.0399933, 189.7200012, 170.3399963, 175.4400024, 175.4400024,
+                             185.6399994, 187.6799927, 168.3000031, 176.4600067, 177.4799957,
+                             185.6399994, 167.2799988, 178.5, 182.5800018, 179.5200043,
+                             173.3999939, 185.6399994, 196.8600006, 183.6000061, 162.1799927,
+                             176.4600067, 189.7200012, 208.0800018, 177.4799957, 163.1999969,
+                             187.6799927, 196.8600006, 180.5399933, 188.6999969, 163.1999969,
+                             187.6799927, 168.3000031, 182.5800018, 181.5599976, 174.4199982,
+                             181.5599976, 161.1600037, 163.1999969, 184.6199951, 190.7400055,
+                             181.5599976, 185.6399994, 186.6600037, 173.3999939, 172.3800049,
+                             179.5200043, 187.6799927, 191.7599945, 190.7400055, 183.6000061,
+                             166.2599945, 196.8600006, 172.3800049, 174.4199982, 181.5599976,
+                             177.4799957, 176.4600067, 188.6999969, 184.6199951, 169.3200073,
+                             178.5, 186.6600037, 174.4199982, 185.6399994, 201.9600067,
+                             171.3600006, 177.4799957, 183.6000061, 165.2400055, 189.7200012,
+                             188.6999969, 178.5, 163.1999969, 169.3200073, 178.5,
+                             182.5800018, 173.3999939, 177.4799957, 165.2400055, 163.1999969,
+                             175.4400024, 184.6199951, 189.7200012, 186.6600037, 188.6999969,
+                             163.1999969, 158.1000061, 172.3800049, 186.6600037, 173.3999939,
+                             157.0800018, 158.1000061, 172.3800049, 197.8800049, 171.3600006,
+                             172.3800049, 184.6199951, 173.3999939, 174.4199982, 175.4400024,
+                             166.2599945, 166.2599945, 172.3800049, 171.3600006, 181.5599976,
+                             181.5599976, 187.6799927, 180.5399933, 169.3200073, 182.5800018,
+                             178.5, 179.5200043, 184.6199951, 175.4400024, 175.4400024,
+                             158.1000061, 182.5800018, 196.8600006, 167.2799988, 178.5,
+                             174.4199982, 180.5399933, 195.8399963, 183.6000061, 200.9400024,
+                             189.7200012, 186.6600037, 173.3999939, 173.3999939, 180.5399933,
+                             172.3800049, 157.0800018, 163.1999969, 171.3600006, 190.7400055,
+                             196.8600006, 179.5200043, 175.4400024, 169.3200073, 158.1000061,
+                             157.0800018, 180.5399933, 173.3999939, 170.3399963, 175.4400024,
+                             193.8000031, 170.3399963, 164.2200012, 174.4199982, 185.6399994,
+                             178.5, 176.4600067, 176.4600067, 179.5200043, 176.4600067,
+                             171.3600006, 205.0200043, 184.6199951, 180.5399933, 165.2400055,
+                             167.2799988, 162.1799927, 165.2400055, 180.5399933, 169.3200073,
+                             176.4600067, 182.5800018, 182.5800018, 175.4400024, 186.6600037,
+                             182.5800018, 183.6000061, 163.1999969, 161.1600037, 189.7200012,
+                             181.5599976, 187.6799927, 173.3999939, 173.3999939, 177.4799957,
+                             179.5200043, 198.8999939, 177.4799957, 183.6000061, 154.0200043,
+                             188.6999969, 181.5599976, 177.4799957, 174.4199982, 202.9799957,
+                             168.3000031, 164.2200012, 187.6799927, 171.3600006, 189.7200012,
+                             185.6399994, 187.6799927, 157.0800018, 193.8000031, 160.1399994,
+                             166.2599945, 193.8000031, 166.2599945, 168.3000031, 179.5200043,
+                             181.5599976, 172.3800049, 183.6000061, 184.6199951, 180.5399933,
+                             177.4799957, 192.7799988, 171.3600006, 197.8800049, 190.7400055,
+                             182.5800018, 178.5, 189.7200012, 172.3800049, 199.9199982,
+                             183.6000061, 179.5200043, 170.3399963, 179.5200043, 181.5599976,
+                             178.5, 186.6600037, 177.4799957, 160.1399994, 176.4600067,
+                             173.3999939, 168.3000031, 180.5399933, 179.5200043, 175.4400024,
+                             188.6999969, 175.4400024, 178.5, 161.1600037, 181.5599976,
+                             184.6199951, 169.3200073, 187.6799927, 164.2200012, 176.4600067,
+                             176.4600067, 174.4199982, 189.7200012, 192.7799988, 181.5599976,
+                             165.2400055, 173.3999939, 184.6199951, 164.2200012, 181.5599976,
+                             167.2799988, 157.0800018, 175.4400024, 172.3800049, 172.3800049,
+                             170.3399963, 166.2599945, 185.6399994, 175.4400024, 184.6199951,
+                             179.5200043, 198.8999939, 189.7200012, 164.2200012, 198.8999939,
+                             169.3200073, 183.6000061, 191.7599945, 168.3000031, 178.5,
+                             172.3800049, 169.3200073, 196.8600006, 170.3399963, 192.7799988,
+                             183.6000061, 186.6600037, 181.5599976, 187.6799927, 198.8999939,
+                             167.2799988, 177.4799957, 165.2400055, 173.3999939, 182.5800018,
+                             190.7400055, 167.2799988, 184.6199951, 180.5399933, 165.2400055,
+                             166.2599945, 162.1799927, 175.4400024, 169.3200073, 187.6799927,
+                             155.0399933, 173.3999939, 165.2400055, 174.4199982, 183.6000061,
+                             167.2799988, 186.6600037, 175.4400024, 173.3999939, 177.4799957,
+                             192.7799988, 180.5399933, 191.7599945, 185.6399994, 194.8200073,
+                             201.9600067, 166.2599945, 171.3600006, 177.4799957, 194.8200073,
+                             191.7599945, 177.4799957, 167.2799988, 188.6999969, 172.3800049,
+                             162.1799927, 169.3200073, 198.8999939, 183.6000061, 170.3399963,
+                             190.7400055, 170.3399963, 169.3200073, 185.6399994, 181.5599976,
+                             166.2599945, 187.6799927, 169.3200073, 157.0800018, 165.2400055,
+                             176.4600067, 174.4199982, 166.2599945, 177.4799957, 195.8399963,
+                             187.6799927, 186.6600037, 194.8200073, 181.5599976, 172.3800049,
+                             166.2599945, 168.3000031, 183.6000061, 168.3000031, 174.4199982,
+                             185.6399994, 180.5399933, 181.5599976, 189.7200012, 172.3800049,
+                             183.6000061, 187.6799927, 183.6000061, 200.9400024, 184.6199951,
+                             173.3999939, 176.4600067, 172.3800049, 169.3200073, 166.2599945,
+                             186.6600037, 181.5599976, 161.1600037, 182.5800018, 179.5200043,
+                             178.5, 174.4199982, 170.3399963, 179.5200043, 193.8000031,
+                             188.6999969, 146.8800049, 192.7799988, 171.3600006, 178.5,
+                             177.4799957, 184.6199951, 180.5399933, 163.1999969, 159.1199951,
+                             160.1399994, 178.5, 176.4600067, 176.4600067, 192.7799988,
+                             161.1600037, 166.2599945, 162.1799927, 172.3800049, 175.4400024,
+                             168.3000031, 201.9600067, 188.6999969, 185.6399994, 175.4400024,
+                             175.4400024, 182.5800018, 182.5800018, 172.3800049, 175.4400024,
+                             179.5200043, 184.6199951, 163.1999969, 195.8399963, 180.5399933,
+                             170.3399963, 212.1600037, 166.2599945, 187.6799927, 179.5200043,
+                             178.5, 176.4600067, 172.3800049, 183.6000061, 179.5200043,
+                             176.4600067, 185.6399994, 161.1600037, 187.6799927, 167.2799988,
+                             187.6799927, 199.9199982, 187.6799927, 169.3200073, 158.1000061,
+                             200.9400024, 191.7599945, 179.5200043, 170.3399963, 186.6600037,
+                             170.3399963, 184.6199951, 189.7200012, 197.8800049, 186.6600037,
+                             171.3600006, 164.2200012, 183.6000061, 180.5399933, 165.2400055,
+                             160.1399994, 183.6000061, 166.2599945, 183.6000061, 196.8600006,
+                             175.4400024, 172.3800049, 172.3800049, 181.5599976, 177.4799957,
+                             173.3999939, 176.4600067, 180.5399933, 176.4600067, 178.5,
+                             163.1999969, 189.7200012, 175.4400024, 174.4199982, 185.6399994,
+                             182.5800018, 169.3200073, 194.8200073, 192.7799988, 188.6999969,
+                             193.8000031, 175.4400024, 165.2400055, 180.5399933, 184.6199951,
+                             176.4600067, 171.3600006, 188.6999969, 199.9199982, 183.6000061,
+                             169.3200073, 175.4400024, 180.5399933, 174.4199982, 167.2799988,
+                             184.6199951, 177.4799957, 180.5399933, 180.5399933, 184.6199951,
+                             178.5, 186.6600037, 161.1600037, 183.6000061, 168.3000031,
+                             188.6999969, 184.6199951, 171.3600006, 185.6399994, 167.2799988,
+                             162.1799927, 186.6600037, 175.4400024, 166.2599945, 182.5800018,
+                             171.3600006, 176.4600067, 192.7799988, 178.5, 168.3000031,
+                             182.5800018, 184.6199951, 148.9199982, 172.3800049, 168.3000031,
+                             181.5599976, 154.0200043, 166.2599945, 174.4199982, 166.2599945,
+                             178.5, 184.6199951, 176.4600067, 171.3600006, 188.6999969,
+                             177.4799957, 172.3800049, 176.4600067, 179.5200043, 176.4600067,
+                             161.1600037, 166.2599945, 181.5599976, 170.3399963, 176.4600067,
+                             163.1999969, 183.6000061, 206.0399933, 171.3600006, 186.6600037,
+                             176.4600067, 163.1999969, 179.5200043, 176.4600067, 182.5800018,
+                             167.2799988, 174.4199982, 188.6999969, 177.4799957, 190.7400055,
+                             176.4600067, 175.4400024, 174.4199982, 179.5200043, 182.5800018,
+                             182.5800018, 175.4400024, 165.2400055, 180.5399933, 178.5,
+                             176.4600067, 191.7599945, 173.3999939, 186.6600037, 165.2400055,
+                             173.3999939, 188.6999969, 172.3800049, 168.3000031, 193.8000031,
+                             180.5399933, 172.3800049, 180.5399933, 189.7200012, 170.3399963,
+                             170.3399963, 195.8399963, 173.3999939, 187.6799927, 165.2400055,
+                             147.8999939, 173.3999939, 184.6199951, 171.3600006, 178.5,
+                             163.1999969, 171.3600006, 190.7400055, 180.5399933, 156.0599976,
+                             164.2200012, 172.3800049, 202.9799957, 201.9600067, 175.4400024,
+                             193.8000031, 176.4600067, 182.5800018, 184.6199951, 178.5,
+                             169.3200073, 166.2599945, 173.3999939, 182.5800018, 171.3600006,
+                             167.2799988, 185.6399994, 179.5200043, 188.6999969, 194.8200073,
+                             177.4799957, 177.4799957, 164.2200012, 189.7200012, 193.8000031,
+                             173.3999939, 178.5, 184.6199951, 191.7599945, 191.7599945,
+                             162.1799927, 171.3600006, 165.2400055, 197.8800049, 169.3200073,
+                             178.5, 174.4199982, 183.6000061, 189.7200012, 181.5599976,
+                             185.6399994, 179.5200043, 181.5599976, 161.1600037, 173.3999939,
+                             187.6799927, 177.4799957, 182.5800018, 171.3600006, 188.6999969,
+                             194.8200073, 176.4600067, 182.5800018, 157.0800018, 169.3200073,
+                             170.3399963, 168.3000031, 190.7400055, 170.3399963, 185.6399994,
+                             195.8399963, 172.3800049, 213.1799927, 187.6799927, 187.6799927,
+                             178.5, 181.5599976, 183.6000061, 163.1999969, 169.3200073,
+                             180.5399933, 204, 175.4400024, 158.1000061, 174.4199982,
+                             183.6000061, 181.5599976, 181.5599976, 180.5399933, 172.3800049,
+                             188.6999969, 166.2599945, 196.8600006, 187.6799927, 196.8600006,
+                             164.2200012, 160.1399994, 165.2400055, 172.3800049, 177.4799957,
+                             173.3999939, 168.3000031, 181.5599976, 164.2200012, 177.4799957,
+                             182.5800018, 162.1799927, 168.3000031, 181.5599976, 179.5200043,
+                             171.3600006, 173.3999939, 176.4600067, 197.8800049, 172.3800049,
+                             189.7200012, 172.3800049, 161.1600037, 183.6000061, 173.3999939,
+                             188.6999969, 161.1600037, 180.5399933, 160.1399994, 166.2599945,
+                             164.2200012, 182.5800018, 169.3200073, 176.4600067, 186.6600037,
+                             192.7799988, 163.1999969, 153, 182.5800018, 184.6199951,
+                             166.2599945, 167.2799988, 186.6600037, 177.4799957, 174.4199982,
+                             175.4400024, 164.2200012, 180.5399933, 178.5, 176.4600067,
+                             164.2200012, 160.1399994, 178.5, 176.4600067, 180.5399933,
+                             182.5800018, 179.5200043, 169.3200073, 178.5, 172.3800049,
+                             171.3600006, 178.5, 166.2599945, 171.3600006, 181.5599976,
+                             162.1799927, 157.0800018, 189.7200012, 177.4799957, 164.2200012,
+                             179.5200043, 187.6799927, 172.3800049, 158.1000061, 164.2200012,
+                             177.4799957, 165.2400055, 185.6399994, 174.4199982, 188.6999969,
+                             173.3999939, 168.3000031, 164.2200012, 187.6799927, 209.1000061,
+                             172.3800049, 177.4799957, 170.3399963, 172.3800049, 164.2200012,
+                             190.7400055, 177.4799957, 165.2400055, 168.3000031, 173.3999939,
+                             0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                             0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                             0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                             0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                             0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                             0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                             0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                             0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                             0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                             0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                             0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                             0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                             0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
+                         };
+
 int main (void)
 {
     plan_tests(26);
 
-    diag("psStats Tests with sample SDSS data from RHL");
+    diag("psStats Tests with sample SDSS data from RHL and Megacam from EAM");
     diag("this file does not yet define a specific test");
-    diag ("the fitted mean is currently wrong for these two data sets");
+    diag("the fitted mean is currently wrong for these two data sets");
 
     {
@@ -243 +515 @@
 
         psVectorStats (stats, y, NULL, NULL, 1);
-        ok (1, "fitted  mean   %f, stdev %f", stats->fittedMean,   stats->fittedStdev);
-        ok (1, "robust  median %f, stdev %f", stats->robustMedian, stats->robustStdev);
-        ok (1, "clipped mean   %f, stdev %f", stats->clippedMean,  stats->clippedStdev);
-        ok (1, "sample  mean   %f, stdev %f", stats->sampleMean,   stats->sampleStdev);
-        ok (1, "sample  median %f, stdev %f", stats->sampleMedian, stats->sampleStdev);
+        ok (1, "sample  mean    %f, stdev %f", stats->sampleMean,   stats->sampleStdev);
+        ok (1, "sample  median  %f, stdev %f", stats->sampleMedian, stats->sampleStdev);
+        ok (1, "clipped mean    %f, stdev %f", stats->clippedMean,  stats->clippedStdev);
+        ok (1, "robust  median  %f, stdev %f", stats->robustMedian, stats->robustStdev);
+        ok (1, "fitted  mean    %f, stdev %f", stats->fittedMean,   stats->fittedStdev);
+        psFree (stats);
+
+        stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2 | PS_STAT_USE_BINSIZE);
+        stats->binsize = 1.0;
+        psVectorStats (stats, y, NULL, NULL, 1);
+        ok (1, "fitted  mean v2 %f, stdev %f", stats->fittedMean,   stats->fittedStdev);
+        psFree (stats);
 
         psFree (y);
-        psFree (stats);
         ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
     }
@@ -273 +551 @@
 
         psVectorStats (stats, y, NULL, NULL, 1);
-        ok (1, "fitted  mean   %f, stdev %f", stats->fittedMean,   stats->fittedStdev);
-        ok (1, "robust  median %f, stdev %f", stats->robustMedian, stats->robustStdev);
-        ok (1, "clipped mean   %f, stdev %f", stats->clippedMean,  stats->clippedStdev);
         ok (1, "sample  mean   %f, stdev %f", stats->sampleMean,   stats->sampleStdev);
         ok (1, "sample  median %f, stdev %f", stats->sampleMedian, stats->sampleStdev);
+        ok (1, "clipped mean   %f, stdev %f", stats->clippedMean,  stats->clippedStdev);
+        ok (1, "robust  median %f, stdev %f", stats->robustMedian, stats->robustStdev);
+        ok (1, "fitted  mean   %f, stdev %f", stats->fittedMean,   stats->fittedStdev);
+        psFree (stats);
+
+        stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2 | PS_STAT_USE_BINSIZE);
+        stats->binsize = 1.0;
+        psVectorStats (stats, y, NULL, NULL, 1);
+        ok (1, "fitted  mean v2 %f, stdev %f", stats->fittedMean,   stats->fittedStdev);
+        psFree (stats);
 
         psFree (y);
-        psFree (stats);
         ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
     }
 
+    {
+        psMemId id = psMemGetId();
+
+        diag("sample 3");
+        psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV |
+                                       PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV |
+                                       PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV |
+                                       PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN |
+                                       PS_STAT_SAMPLE_STDEV | PS_STAT_USE_BINSIZE);
+        stats->binsize = 1.0;
+
+
+        // copy data in static array
+        int nPts = sizeof(yraw_03) / sizeof (float);
+        psVector *y = psVectorAlloc (nPts, PS_TYPE_F32);
+        for (int i = 0; i < y->n; i++) {
+            y->data.F32[i] = yraw_03[i];
+        }
+
+        psVectorStats (stats, y, NULL, NULL, 1);
+        ok (1, "sample  mean   %f, stdev %f", stats->sampleMean,   stats->sampleStdev);
+        ok (1, "sample  median %f, stdev %f", stats->sampleMedian, stats->sampleStdev);
+        ok (1, "clipped mean   %f, stdev %f", stats->clippedMean,  stats->clippedStdev);
+        ok (1, "robust  median %f, stdev %f", stats->robustMedian, stats->robustStdev);
+        ok (1, "fitted  mean   %f, stdev %f", stats->fittedMean,   stats->fittedStdev);
+        psFree (stats);
+
+        stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2 | PS_STAT_USE_BINSIZE);
+        stats->binsize = 1.0;
+        psVectorStats (stats, y, NULL, NULL, 1);
+        ok (1, "fitted  mean v2 %f, stdev %f", stats->fittedMean,   stats->fittedStdev);
+        psFree (stats);
+
+        psFree (y);
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
     return exit_status();
 }