Index: trunk/psLib/src/math/psStats.c =================================================================== --- trunk/psLib/src/math/psStats.c (revision 36238) +++ trunk/psLib/src/math/psStats.c (revision 36290) @@ -1049,24 +1049,24 @@ #if (CZW) - printf("CZW: bad sigma?: %f %f %f %f %f %f %f %f %f %f\n", - binH2F32,binL2F32,binHiF32,binLoF32,binH4F32,binL4F32, - sigma1,sigma2,sigma4,sigma); + // printf("CZW: bad sigma?: %f %f %f %f %f %f %f %f %f %f\n", + // binH2F32,binL2F32,binHiF32,binLoF32,binH4F32,binL4F32, + // sigma1,sigma2,sigma4,sigma); - printf("CZW (%d): median %f sigma %f delta: %f \n\t %f %f %f %f %f %f %f \n\t %f %f %f %f %f %f %f\n", + printf("CZW Robust (%d): median %f sigma %f delta: %f \n\t %f %f %f %f %f %f %f \n\t %f %f %f %f %f %f %f\n", iterate, - stats->robustMedian,stats->robustStdev, - fabs(cumulative->bounds->data.F32[binMedian] - cumulative->bounds->data.F32[binMedian + 1]), - - cumulative->bounds->data.F32[binMedian-3],cumulative->bounds->data.F32[binMedian-2], - cumulative->bounds->data.F32[binMedian-1], - cumulative->bounds->data.F32[binMedian], - cumulative->bounds->data.F32[binMedian+1], - cumulative->bounds->data.F32[binMedian+2],cumulative->bounds->data.F32[binMedian+3], - - cumulative->nums->data.F32[binMedian-3],cumulative->nums->data.F32[binMedian-2], - cumulative->nums->data.F32[binMedian-1], - cumulative->nums->data.F32[binMedian], - cumulative->nums->data.F32[binMedian+1], - cumulative->nums->data.F32[binMedian+2],cumulative->nums->data.F32[binMedian+3]); + stats->robustMedian,stats->robustStdev, + fabs(cumulative->bounds->data.F32[binMedian] - cumulative->bounds->data.F32[binMedian + 1]), + + cumulative->bounds->data.F32[binMedian-3],cumulative->bounds->data.F32[binMedian-2], + cumulative->bounds->data.F32[binMedian-1], + cumulative->bounds->data.F32[binMedian], + cumulative->bounds->data.F32[binMedian+1], + cumulative->bounds->data.F32[binMedian+2],cumulative->bounds->data.F32[binMedian+3], + + cumulative->nums->data.F32[binMedian-3],cumulative->nums->data.F32[binMedian-2], + cumulative->nums->data.F32[binMedian-1], + cumulative->nums->data.F32[binMedian], + cumulative->nums->data.F32[binMedian+1], + cumulative->nums->data.F32[binMedian+2],cumulative->nums->data.F32[binMedian+3]); // PS_VECTOR_PRINT_F32(histogram->bounds); // PS_VECTOR_PRINT_F32(histogram->nums); @@ -1262,6 +1262,15 @@ // 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; + // psF32 dN = PS_MAX (1, PS_MIN (4, stats->robustN50 / 8)); + + // CZW 2013-11-20: We know that the histogram is going to be basically Gaussian. + // Furthermore, we only use the inner +/- 2 sigma parts. Therefore, define the + // binsize such that the bin at 2 sigma contains ~50 points (S/N ~ 7). robustN50 + // contains half the total points, so 2 * robustN50 / 50 is the fraction of all + // points in the 2 sigma bin. Dance the erf() relations around, and it looks like + // there's a factor of about 1/20 to include. Keep the PS_MAX to ensure we never bin + // wider than 1 sigma when the number of points is small. + psF32 dN = PS_MAX(1, (stats->robustN50 / 500.0)); + binSize = guessStdev / dN; } @@ -1291,9 +1300,9 @@ // We no longer want to reset the binSize here, as it can cause odd things. Better to select // a number of bins, and then set the min/max values to put those bins sanely around the mean. - long numBins = PS_MAX (50, PS_MIN (10000, (max - min) / binSize)); + // long numBins = PS_MAX (50, PS_MIN (10000, (max - min) / binSize)); // binSize = (max - min) / (float) numBins; 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); + // psTrace(TRACE, 6, "The numBins is %ld\n", numBins); @@ -1305,5 +1314,6 @@ max = guessMean + FITTED_CLIPPING_NUM * guessStdev; } - + long numBins = PS_MAX (50, PS_MIN (10000, (max - min) / binSize)); + if (CZW) { printf("I've clipped: %f %f => %f %f ; %f %ld\n",guessMean,guessStdev,min,max,binSize,stats->robustN50); } psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers) if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) { @@ -1357,7 +1367,16 @@ } psTrace (TRACE, 6, "(%f = %.0f) ", histogram->bounds->data.F32[i], histogram->nums->data.F32[i]); + if (CZW) { printf("CENTERED_HISTOGRAM: %f %f\n", + PS_BIN_MIDPOINT(histogram,i), + histogram->nums->data.F32[i]); } } psTrace (TRACE, 6, "\n"); + if (CZW) { printf("Bin selection done: %ld %f %f %ld %f %f %ld %f %f\n", + binMin,PS_BIN_MIDPOINT(histogram,binMin),histogram->nums->data.F32[binMin], + binMax,PS_BIN_MIDPOINT(histogram,binMax),histogram->nums->data.F32[binMax], + binPeak,PS_BIN_MIDPOINT(histogram,binPeak),histogram->nums->data.F32[binPeak]); + } + // 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); @@ -1411,4 +1430,13 @@ psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); bool status = psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x); +#if (CZW && 1) + printf("CZW: LowfitPoly: %f %f %f\n",poly->coeff[0],poly->coeff[1],poly->coeff[2]); + for (long i = 0; i < x->n; i++) { + printf("CZW: Lowfit: %d %ld %f %f %f\n", + status,i,x->data.F32[i],y->data.F32[i], + poly->coeff[0] + poly->coeff[1] * x->data.F32[i] + + poly->coeff[2] * pow(x->data.F32[i],2)); + } +#endif psFree(x); psFree(y); @@ -1417,5 +1445,4 @@ psErrorClear(); COUNT_WARNING(10, 100, "Failed to fit a gaussian to the robust histogram.\n"); - psFree(poly); psFree(histogram); @@ -1426,4 +1453,5 @@ if (poly->coeff[2] >= 0.0) { COUNT_WARNING(10, 100, "Failed parabolic fit: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]); + psFree(poly); psFree(histogram); @@ -1478,4 +1506,5 @@ } } + psTrace(TRACE, 6, "Lower bound for symmetric range: %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binS), binS); @@ -1499,8 +1528,9 @@ psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); bool status = psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x); -#if (CZW && 0) +#if (CZW && 1) + printf("CZW: FullfitPoly: %f %f %f\n",poly->coeff[0],poly->coeff[1],poly->coeff[2]); for (long i = 0; i < x->n; i++) { - printf("CZW: Dcheck: %ld %f %f %f\n", - i,x->data.F32[i],y->data.F32[i], + printf("CZW: Fullfit: %d %ld %f %f %f\n", + status,i,x->data.F32[i],y->data.F32[i], poly->coeff[0] + poly->coeff[1] * x->data.F32[i] + poly->coeff[2] * pow(x->data.F32[i],2));