Index: /branches/eam_branch_20070830/psLib/src/math/psStats.c =================================================================== --- /branches/eam_branch_20070830/psLib/src/math/psStats.c (revision 14769) +++ /branches/eam_branch_20070830/psLib/src/math/psStats.c (revision 14770) @@ -13,6 +13,6 @@ * use ->min and ->max (PS_STAT_USE_RANGE) * - * @version $Revision: 1.215.2.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-09-02 02:03:58 $ + * @version $Revision: 1.215.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-07 20:16:59 $ * * Copyright 2006 IfA, University of Hawaii @@ -1469,5 +1469,7 @@ valPeak = histogram->nums->data.F32[binPeak]; } - } + psTrace (TRACE, 6, "(%f = %.0f) ", histogram->bounds->data.F32[i], histogram->nums->data.F32[i]); + } + psTrace (TRACE, 6, "\n"); // assume a reasonably well-defined gaussian-like population; run from peak out until val < 0.25*peak @@ -1477,4 +1479,5 @@ psTrace(TRACE, 6, "The clipped max is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMax - 1), binMax - 1); psTrace(TRACE, 6, "The clipped peak is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binPeak), binPeak); + psTrace(TRACE, 6, "The clipped peak value is %f\n", histogram->nums->data.F32[binPeak]); { @@ -1529,5 +1532,5 @@ if (iteration == 0) { guessStdev = 0.25*guessStdev; - psTrace(TRACE, 6, "*** retry stdev is %f.\n", guessStdev); + psTrace(TRACE, 6, "*** retry, new stdev is %f.\n", guessStdev); continue; } @@ -1551,4 +1554,5 @@ } + // for test, measure the same result for the upper section { // fit the upper half of the distribution @@ -1592,5 +1596,5 @@ } - // calculate upper mean & stdev from parabolic fit -- use this as the result + // calculate upper mean & stdev from parabolic fit -- ignore this value #ifndef PS_NO_TRACE float upperStdev = sqrt(-0.5/poly->coeff[2]); @@ -1601,4 +1605,10 @@ #endif + // if the resulting value is outside of the range binMin - binMax, use the upper value + if (done && (guessMean > PS_BIN_MIDPOINT(histogram, binMax - 1))) { + guessMean = upperMean; + guessStdev = upperStdev; + } + psFree (poly); } @@ -1619,4 +1629,296 @@ stats->results |= PS_STAT_FITTED_MEAN_V3; stats->results |= PS_STAT_FITTED_STDEV_V3; + + return true; +} + +/******************** + * perform an asymmetric fit to the population + * vectorFittedStats_v4 requires guess for fittedMean and fittedStdev + * robustN50 should also be set + * gaussian fit is performed using 2D polynomial to ln(y) + * this version follows the upper portion of the distribution until it passes 0.5*peak + ********************/ +static bool vectorFittedStats_v4 (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(TRACE, 4, "---- %s() end ----\n", __func__); + return false; + } + + float guessStdev = stats->robustStdev; // pass the guess sigma + float guessMean = stats->robustMedian; // pass the guess mean + + psTrace(TRACE, 6, "The ** starting ** guess mean is %f.\n", guessMean); + psTrace(TRACE, 6, "The ** starting ** 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 + + psF32 binSize = 1; + if (stats->options & PS_STAT_USE_BINSIZE) { + // Set initial bin size to the specified value. + binSize = stats->binsize; + psTrace(TRACE, 6, "Setting initial robust bin size to %.2f\n", binSize); + } else { + // construct a histogram with (sigma/2 < binsize < sigma) + // set roughly so that the lowest bins have about 2 cnts + // Nsmallest ~ N50 / (4*dN)) + psF32 dN = PS_MAX (1, PS_MIN (4, 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(TRACE, 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(TRACE, 6, "The new min/max values are (%f, %f).\n", min, max); + psTrace(TRACE, 6, "The new bin size is %f.\n", binSize); + psTrace(TRACE, 6, "The numBins is %ld\n", numBins); + + psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers) + if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) { + psError(PS_ERR_UNKNOWN, false, "Unable to generate histogram for fitted statistics v4.\n"); + psFree(histogram); + psFree(statsMinMax); + return false; + } + if (psTraceGetLevel("psLib.math") >= 8) { + PS_VECTOR_PRINT_F32(histogram->nums); + } + + // now fit a Gaussian to the upper and lower halves about the peak independently + + // set the full-range upper and lower limits + 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); + } + + // select the min and max bins, saturating on the lower and upper end-points + long binMin, binMax; + PS_BIN_FOR_VALUE (binMin, histogram->bounds, guessMean - minFitSigma*guessStdev, 0); + PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, 0); + if (binMin == binMax) { + psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n"); + psFree(statsMinMax); + psTrace(TRACE, 4, "---- %s(false) end ----\n", __func__); + return false; + } + + // search for mode (peak of histogram within range mean-2sigma - mean+2sigma + long binPeak = binMin; + float valPeak = histogram->nums->data.F32[binPeak]; + for (int i = binMin; i < binMax; i++) { + if (histogram->nums->data.F32[i] > valPeak) { + binPeak = i; + valPeak = histogram->nums->data.F32[binPeak]; + } + psTrace (TRACE, 6, "(%f = %.0f) ", histogram->bounds->data.F32[i], histogram->nums->data.F32[i]); + } + psTrace (TRACE, 6, "\n"); + + // assume a reasonably well-defined gaussian-like population; run from peak out until val < 0.25*peak + + psTrace(TRACE, 6, "The clipped numBins is %ld\n", binMax - binMin); + psTrace(TRACE, 6, "The clipped min is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMin), binMin); + psTrace(TRACE, 6, "The clipped max is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMax - 1), binMax - 1); + psTrace(TRACE, 6, "The clipped peak is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binPeak), binPeak); + psTrace(TRACE, 6, "The clipped peak value is %f\n", histogram->nums->data.F32[binPeak]); + + { + // fit the lower half of the distribution + // run down until we drop below 0.25*valPeak + // run up until we drop below 0.50*valPeak + long binS = binMin; + long binE = PS_MIN (binPeak + 3, binMax); + for (int i = binPeak - 3; i >= binMin; i--) { + if (histogram->nums->data.F32[i] < 0.25*valPeak) { + binS = i; + break; + } + } + for (int i = binPeak + 3; i < binMax; i++) { + if (histogram->nums->data.F32[i] < 0.50*valPeak) { + binE = i; + break; + } + } + psTrace(TRACE, 6, "Lower bound for lower half: %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binS), binS); + psTrace(TRACE, 6, "Upper bound for lower half: %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binE), binE); + + psVector *y = psVectorAllocEmpty(binE - binS, PS_TYPE_F32); // Vector of coordinates + psVector *x = psVectorAllocEmpty(binE - binS, PS_TYPE_F32); // Vector of ordinates + long j = 0; + for (long i = binS; i < binE; 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]); // note this is the natural log: expected distribution is A exp(-(x-xo)^2/2sigma^2) + j++; + } + y->n = x->n = j; + + // fit 2nd order polynomial to ln(y) = -(x-xo)^2/2sigma^2 + psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); + bool status = psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x); + psFree(x); + psFree(y); + + if (!status) { + psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n"); + psFree(poly); + psFree(histogram); + psFree(statsMinMax); + psTrace(TRACE, 4, "---- %s(false) end ----\n", __func__); + return false; + } + + if (poly->coeff[2] >= 0.0) { + psTrace(TRACE, 6, "Failed parabolic fit: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]); + psFree(poly); + psFree(histogram); + psFree(statsMinMax); + + // sometimes, the guessStdev is much too large. in this case, the entire real population tends to be found + // in a single bin. make one attempt to recover by dropping the guessStdev down by a jump and trying again + if (iteration == 0) { + guessStdev = 0.25*guessStdev; + psTrace(TRACE, 6, "*** retry, new stdev is %f.\n", guessStdev); + continue; + } + + psError(PS_ERR_UNKNOWN, false, "fit did not converge\n"); + psTrace(TRACE, 4, "---- %s(false) end ----\n", __func__); + return false; + } + + // calculate lower mean & stdev from parabolic fit -- use this as the result + guessStdev = sqrt(-0.5/poly->coeff[2]); + guessMean = poly->coeff[1]*PS_SQR(guessStdev); + if (guessStdev > 0.75*stats->robustStdev) { + done = true; + } + psTrace(TRACE, 6, "Parabolic Lower fit results: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]); + psTrace(TRACE, 6, "The lower mean is %f.\n", guessMean); + psTrace(TRACE, 6, "The lower stdev is %f.\n", guessStdev); + + psFree(poly); + } + + // if we converge on a solution outside the range binMin - binMax, use a more conservative range + if (done && (guessMean > PS_BIN_MIDPOINT(histogram, binMax - 1))) { + psTrace(TRACE, 6, "Inconsistent result, re-trying the fit\n"); + + // fit a symmetric distribution + // run up until we drop below 0.15*valPeak + // run up until we drop below 0.15*valPeak + long binS = binMin; + long binE = binMax; + for (int i = binPeak - 3; i >= binMin; i--) { + if (histogram->nums->data.F32[i] < 0.15*valPeak) { + binS = i; + break; + } + } + for (int i = binPeak + 3; i < binMax; i++) { + if (histogram->nums->data.F32[i] < 0.15*valPeak) { + binE = i; + break; + } + } + psTrace(TRACE, 6, "Lower bound for symmetric range: %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binS), binS); + psTrace(TRACE, 6, "Upper bound for symmetric range: %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binE), binE); + + psVector *y = psVectorAllocEmpty(binE - binS, PS_TYPE_F32); // Vector of coordinates + psVector *x = psVectorAllocEmpty(binE - binS, PS_TYPE_F32); // Vector of ordinates + long j = 0; + for (long i = binS; i < binE; 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]); // note this is the natural log: expected distribution is A exp(-(x-xo)^2/2sigma^2) + j++; + } + y->n = x->n = j; + + // fit 2nd order polynomial to ln(y) = -(x-xo)^2/2sigma^2 + psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); + bool status = psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x); + psFree(x); + psFree(y); + + if (!status) { + psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n"); + psFree(poly); + psFree(histogram); + psFree(statsMinMax); + psTrace(TRACE, 4, "---- %s(false) end ----\n", __func__); + return false; + } + + // calculate upper mean & stdev from parabolic fit -- ignore this value +#ifndef PS_NO_TRACE + guessStdev = sqrt(-0.5/poly->coeff[2]); + guessMean = poly->coeff[1]*PS_SQR(guessStdev); + psTrace(TRACE, 6, "Parabolic Symmetric fit results: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]); + psTrace(TRACE, 6, "The symmetric mean is %f.\n", guessMean); + psTrace(TRACE, 6, "The symmetric stdev is %f.\n", guessStdev); +#endif + + psFree (poly); + } + + // Clean up after fitting + psFree (histogram); + psFree (statsMinMax); + } + + // The fitted mean is the Gaussian mean. + stats->fittedMean = guessMean; + psTrace(TRACE, 6, "The fitted mean is %f.\n", stats->fittedMean); + + // The fitted standard deviation + stats->fittedStdev = guessStdev; + psTrace(TRACE, 6, "The fitted stdev is %f.\n", stats->fittedStdev); + + stats->results |= PS_STAT_FITTED_MEAN_V4; + stats->results |= PS_STAT_FITTED_STDEV_V4; return true; @@ -1946,4 +2248,15 @@ // ************************************************************************ + if (stats->options & (PS_STAT_FITTED_MEAN_V4 | PS_STAT_FITTED_STDEV_V4)) { + if (stats->options & (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV)) { + psAbort("you may not specify both FITTED_MEAN and FITTED_MEAN_V4"); + } + if (!vectorFittedStats_v4(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)) { @@ -1992,4 +2305,7 @@ READ_STAT("FITTED_MEAN_V3", PS_STAT_FITTED_MEAN_V3); READ_STAT("FITTED_STDEV_V3", PS_STAT_FITTED_STDEV_V3); + READ_STAT("FITTED_V4", PS_STAT_FITTED_MEAN_V4); + READ_STAT("FITTED_MEAN_V4", PS_STAT_FITTED_MEAN_V4); + READ_STAT("FITTED_STDEV_V4", PS_STAT_FITTED_STDEV_V4); READ_STAT("CLIPPED", PS_STAT_CLIPPED_MEAN); READ_STAT("CLIPPED_MEAN", PS_STAT_CLIPPED_MEAN); @@ -2025,4 +2341,6 @@ WRITE_STAT("FITTED_MEAN_V3", PS_STAT_FITTED_MEAN_V3); WRITE_STAT("FITTED_STDEV_V3", PS_STAT_FITTED_STDEV_V3); + WRITE_STAT("FITTED_MEAN_V4", PS_STAT_FITTED_MEAN_V4); + WRITE_STAT("FITTED_STDEV_V4", PS_STAT_FITTED_STDEV_V4); WRITE_STAT("CLIPPED_MEAN", PS_STAT_CLIPPED_MEAN); WRITE_STAT("CLIPPED_STDEV", PS_STAT_CLIPPED_STDEV); @@ -2079,4 +2397,6 @@ case PS_STAT_FITTED_MEAN_V3: case PS_STAT_FITTED_STDEV_V3: + case PS_STAT_FITTED_MEAN_V4: + case PS_STAT_FITTED_STDEV_V4: case PS_STAT_CLIPPED_MEAN: case PS_STAT_CLIPPED_STDEV: @@ -2121,4 +2441,8 @@ case PS_STAT_FITTED_STDEV_V3: return stats->fittedStdev; + case PS_STAT_FITTED_MEAN_V4: + return stats->fittedMean; + case PS_STAT_FITTED_STDEV_V4: + return stats->fittedStdev; case PS_STAT_CLIPPED_MEAN: return stats->clippedMean; Index: /branches/eam_branch_20070830/psLib/src/math/psStats.h =================================================================== --- /branches/eam_branch_20070830/psLib/src/math/psStats.h (revision 14769) +++ /branches/eam_branch_20070830/psLib/src/math/psStats.h (revision 14770) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.62.2.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-09-02 02:03:58 $ + * @version $Revision: 1.62.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-07 20:16:59 $ * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii */ @@ -51,6 +51,6 @@ PS_STAT_USE_RANGE = 0x100000, ///< Range PS_STAT_USE_BINSIZE = 0x200000, ///< Binsize - PS_STAT_SPARE2 = 0x400000, ///< Spare 2 - PS_STAT_SPARE3 = 0x800000, ///< Spare 3 + PS_STAT_FITTED_MEAN_V4 = 0x400000, ///< Fitted Mean + PS_STAT_FITTED_STDEV_V4 = 0x800000, ///< Fitted Standard Deviation } psStatsOptions;