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