Index: trunk/psLib/src/math/psStats.c =================================================================== --- trunk/psLib/src/math/psStats.c (revision 1057) +++ trunk/psLib/src/math/psStats.c (revision 1071) @@ -847,4 +847,205 @@ /****************************************************************************** +p_psVectorSampleStdev(myVector, maskVector, maskVal, stats): calculates the +stdev of the input vector. +Inputs + myVector + maskVector + maskVal + stats +Returns + NULL + +NOTE: the mean is always calculated exactly. Robust means are never +calculated in this routine. + *****************************************************************************/ +void p_psVectorSampleStdev(const psVector *restrict myVector, + const psVector *restrict maskVector, + unsigned int maskVal, + psStats *stats) +{ + int i = 0; // Loop index variable + int countInt = 0; // # of data points being used + float countFloat = 0.0; // # of data points being used + float mean = 0.0; // The mean + float diff = 0.0; // Used in calculating stdev + float sumSquares = 0.0; // temporary variable + float sumDiffs = 0.0; // temporary variable + float rangeMin = 0.0; // Exclude data below this + float rangeMax = 0.0; // Exclude date above this + + // This procedure requires the mean. If it has not been already + // calculated, then call p_psVectorSampleMean() + if (0 != isnan(stats->sampleMean)) { + p_psVectorSampleMean(myVector, maskVector, maskVal, stats); + } + mean = stats->sampleMean; + + if (stats->options & PS_STAT_USE_RANGE) { + if (maskVector != NULL) { + for (i=0;in;i++) { + if (!(maskVal & maskVector->data.U8[i]) && + (rangeMin <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= rangeMax)) { + diff = myVector->data.F32[i] - mean; + sumSquares+= (diff * diff); + sumDiffs+= diff; + countInt++; + } + } + } else { + for (i=0;in;i++) { + if ((rangeMin <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= rangeMax)) { + diff = myVector->data.F32[i] - mean; + sumSquares+= (diff * diff); + sumDiffs+= diff; + countInt++; + } + } + countInt = myVector->n; + } + } else { + if (maskVector != NULL) { + for (i=0;in;i++) { + if (!(maskVal & maskVector->data.U8[i])) { + diff = myVector->data.F32[i] - mean; + sumSquares+= (diff * diff); + sumDiffs+= diff; + countInt++; + } + } + } else { + for (i=0;in;i++) { + diff = myVector->data.F32[i] - mean; + sumSquares+= (diff * diff); + sumDiffs+= diff; + countInt++; + } + countInt = myVector->n; + } + } + countFloat = (float) countInt; + + #ifdef DARWIN + + stats->sampleStdev = (float) sqrt( (sumSquares-(sumDiffs * + sumDiffs/countFloat))/ (countFloat-1)); + #else + + stats->sampleStdev = sqrtf( (sumSquares-(sumDiffs * + sumDiffs/countFloat))/ (countFloat-1)); + #endif +} + +/****************************************************************************** +p_psVectorClippedStats(myVector, maskVector, maskVal, stats): calculates the +clipped stats (mean or stdev) of the input vector. + +Inputs + myVector + maskVector + maskVal + stats +Returns + NULL + *****************************************************************************/ +void p_psVectorClippedStats(const psVector *restrict myVector, + const psVector *restrict maskVector, + unsigned int maskVal, + psStats *stats) +{ + int i = 0; // Loop index variable + int j = 0; // Loop index variable + float clippedMean = 0.0; // self-explanatory + float clippedStdev = 0.0; // self-explanatory + float oldStanMean = 0.0; // Temporary variable + float oldStanStdev = 0.0; // Temporary variable + psVector *tmpMask = NULL; // Temporary vector + + // Endure that stats->clipIter is within the proper range. + if (!((CLIPPED_NUM_ITER_LB <= stats->clipIter ) && + (stats->clipIter <= CLIPPED_NUM_ITER_UB))) { + psAbort(__func__, "Unallowed value for clipIter (%d).\n", + stats->clipIter); + } + + // Endure that stats->clipSigma is within the proper range. + if (!((CLIPPED_SIGMA_LB <= stats->clipSigma ) && + (stats->clipSigma <= CLIPPED_SIGMA_UB))) { + psAbort(__func__, "Unallowed value for clipSigma (%f).\n", + stats->clipSigma); + } + + // We allocate a temporary mask vector since during the iterative + // steps that follow, we will be masking off additional data points. + // However, we do no want to modify the original mask vector. + tmpMask = psVectorAlloc(myVector->n, PS_TYPE_U8); + tmpMask->n = myVector->n; + + // If we were called with a mask vector, then initialize the temporary + // mask vector with those values. + if (maskVector != NULL) { + for (i=0;in;i++) { + tmpMask->data.U8[i] = maskVector->data.U8[i]; + } + } + + // 1. Compute the sample median. + // NOTE: This seems odd. Verify with IfA that we want to calculate the + // median here, not the mean. + p_psVectorSampleMedian(myVector, maskVector, maskVal, stats); + + // 2. Compute the sample standard deviation. + p_psVectorSampleStdev(myVector, maskVector, maskVal, stats); + + // 3. Use the sample median as the first estimator of the mean X. + clippedMean = stats->sampleMean; + + // 4. Use the sample stdev as the first estimator of the mean stdev. + clippedStdev = stats->sampleStdev; + + // Must save the old sampleMean and sampleStdev since the following code + // block overwrites them. + oldStanMean = stats->sampleMean; + oldStanStdev = stats->sampleStdev; + + // 5. Repeat N times: + for (i=0;iclipIter;i++) { + for (j=0;jn;j++) { + // a) Exclude all values x_i for which |x_i - x| > K * stdev + if (fabs(myVector->data.F32[j] - clippedMean) > + (stats->clipSigma * clippedStdev)) { + tmpMask->data.U8[i] = 0xff; + } + // b) compute new mean and stdev + p_psVectorSampleMedian(myVector, tmpMask, maskVal, stats); + p_psVectorSampleStdev(myVector, tmpMask, maskVal, stats); + + // c) Use the new mean for x + clippedMean = stats->sampleMean; + + // d) Use the new stdev for stdev + clippedStdev = stats->sampleStdev; + } + } + stats->sampleMean = oldStanMean; + stats->sampleStdev= oldStanStdev; + + // 7. The last calcuated value of x is the cliped mean. + if (stats->options & PS_STAT_CLIPPED_MEAN) { + stats->clippedMean = clippedMean; + } + + // 8. The last calcuated value of stdev is the cliped stdev. + if (stats->options & PS_STAT_CLIPPED_STDEV) { + stats->clippedStdev = clippedStdev; + } + + psVectorFree(tmpMask); +} + + +/****************************************************************************** p_psVectorRobustStats(myVector, maskVector, maskVal, stats): this procedure calculates a variety of robust stat measures: @@ -884,4 +1085,5 @@ float dL = 0.0; int numBins = 0; + psStats *tmpStats = psStatsAlloc(PS_STAT_CLIPPED_STDEV|PS_STAT_CLIPPED_MEAN); // NOTE: The SDRS states that the sample quartiles must be used to @@ -890,17 +1092,41 @@ // size of the data set. We should consult with IfA to ensure that this // is really required. - - if (isnan(stats->sampleUQ) || + /* + if (isnan(stats->sampleUQ) || isnan(stats->sampleLQ)) { - stats->options = stats->options | PS_STAT_SAMPLE_QUARTILE; - p_psVectorSampleQuartiles(myVector, - maskVector, - maskVal, - stats); - } - + stats->options = stats->options | PS_STAT_SAMPLE_QUARTILE; + p_psVectorSampleQuartiles(myVector, + maskVector, + maskVal, + stats); + } + */ // Compute the initial bin size of the robust histogram. - sigmaE = (stats->sampleUQ - stats->sampleLQ) / 1.34f; - binSize = sigmaE / 10.0f; + p_psVectorClippedStats(myVector, maskVector, maskVal, tmpStats); + binSize = stats->clippedStdev / 10.0f; + + // If stats->clippedStdev == 0.0, then all data elements have the same + // value. Therefore, we can set the appropiate results and return. + if (fabs(binSize) <= FLT_EPSILON) { + if (stats->options & PS_STAT_ROBUST_MEAN) { + stats->robustMean = stats->clippedMean; + } + if (stats->options & PS_STAT_ROBUST_MEDIAN) { + stats->robustMedian = stats->clippedMean; + } + if (stats->options & PS_STAT_ROBUST_MODE) { + stats->robustMode = stats->clippedMean; + } + if (stats->options & PS_STAT_ROBUST_STDEV) { + stats->robustStdev = 0.0; + } + if (stats->options & PS_STAT_ROBUST_QUARTILE) { + stats->robustUQ = stats->clippedMean; + stats->robustLQ = stats->clippedMean; + } + psStatsFree(tmpStats); + psHistogramFree(robustHistogram); + return; + } // Detemine minimum and maximum values in the data vector. @@ -912,5 +1138,6 @@ } - // Create the histogram structure (yes, 2 is necessary, not 1). + // Create the histogram structure (yes, 2 is necessary, not 1). Also, + // if we get here, we know that binSize != 0.0. numBins = 2 + (int) ((stats->max - stats->min) / binSize); @@ -973,205 +1200,6 @@ stats->robustNfit = 0.0; stats->robustN50 = 0.0; + psStatsFree(tmpStats); psHistogramFree(robustHistogram); -} - -/****************************************************************************** -p_psVectorSampleStdev(myVector, maskVector, maskVal, stats): calculates the -stdev of the input vector. -Inputs - myVector - maskVector - maskVal - stats -Returns - NULL - -NOTE: the mean is always calculated exactly. Robust means are never -calculated in this routine. - *****************************************************************************/ -void p_psVectorSampleStdev(const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats) -{ - int i = 0; // Loop index variable - int countInt = 0; // # of data points being used - float countFloat = 0.0; // # of data points being used - float mean = 0.0; // The mean - float diff = 0.0; // Used in calculating stdev - float sumSquares = 0.0; // temporary variable - float sumDiffs = 0.0; // temporary variable - float rangeMin = 0.0; // Exclude data below this - float rangeMax = 0.0; // Exclude date above this - - // This procedure requires the mean. If it has not been already - // calculated, then call p_psVectorSampleMean() - if (0 != isnan(stats->sampleMean)) { - p_psVectorSampleMean(myVector, maskVector, maskVal, stats); - } - mean = stats->sampleMean; - - if (stats->options & PS_STAT_USE_RANGE) { - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i]) && - (rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - diff = myVector->data.F32[i] - mean; - sumSquares+= (diff * diff); - sumDiffs+= diff; - countInt++; - } - } - } else { - for (i=0;in;i++) { - if ((rangeMin <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= rangeMax)) { - diff = myVector->data.F32[i] - mean; - sumSquares+= (diff * diff); - sumDiffs+= diff; - countInt++; - } - } - countInt = myVector->n; - } - } else { - if (maskVector != NULL) { - for (i=0;in;i++) { - if (!(maskVal & maskVector->data.U8[i])) { - diff = myVector->data.F32[i] - mean; - sumSquares+= (diff * diff); - sumDiffs+= diff; - countInt++; - } - } - } else { - for (i=0;in;i++) { - diff = myVector->data.F32[i] - mean; - sumSquares+= (diff * diff); - sumDiffs+= diff; - countInt++; - } - countInt = myVector->n; - } - } - countFloat = (float) countInt; - - #ifdef DARWIN - - stats->sampleStdev = (float) sqrt( (sumSquares-(sumDiffs * - sumDiffs/countFloat))/ (countFloat-1)); - #else - - stats->sampleStdev = sqrtf( (sumSquares-(sumDiffs * - sumDiffs/countFloat))/ (countFloat-1)); - #endif -} - -/****************************************************************************** -p_psVectorClippedStats(myVector, maskVector, maskVal, stats): calculates the -clipped stats (mean or stdev) of the input vector. - -Inputs - myVector - maskVector - maskVal - stats -Returns - NULL - *****************************************************************************/ -void p_psVectorClippedStats(const psVector *restrict myVector, - const psVector *restrict maskVector, - unsigned int maskVal, - psStats *stats) -{ - int i = 0; // Loop index variable - int j = 0; // Loop index variable - float clippedMean = 0.0; // self-explanatory - float clippedStdev = 0.0; // self-explanatory - float oldStanMean = 0.0; // Temporary variable - float oldStanStdev = 0.0; // Temporary variable - psVector *tmpMask = NULL; // Temporary vector - - // Endure that stats->clipIter is within the proper range. - if (!((CLIPPED_NUM_ITER_LB <= stats->clipIter ) && - (stats->clipIter <= CLIPPED_NUM_ITER_UB))) { - psAbort(__func__, "Unallowed value for clipIter (%d).\n", - stats->clipIter); - } - - // Endure that stats->clipSigma is within the proper range. - if (!((CLIPPED_SIGMA_LB <= stats->clipSigma ) && - (stats->clipSigma <= CLIPPED_SIGMA_UB))) { - psAbort(__func__, "Unallowed value for clipSigma (%f).\n", - stats->clipSigma); - } - - // We allocate a temporary mask vector since during the iterative - // steps that follow, we will be masking off additional data points. - // However, we do no want to modify the original mask vector. - tmpMask = psVectorAlloc(myVector->n, PS_TYPE_U8); - tmpMask->n = myVector->n; - - // If we were called with a mask vector, then initialize the temporary - // mask vector with those values. - if (maskVector != NULL) { - for (i=0;in;i++) { - tmpMask->data.U8[i] = maskVector->data.U8[i]; - } - } - - // 1. Compute the sample median. - // NOTE: This seems odd. Verify with IfA that we want to calculate the - // median here, not the mean. - p_psVectorSampleMedian(myVector, maskVector, maskVal, stats); - - // 2. Compute the sample standard deviation. - p_psVectorSampleStdev(myVector, maskVector, maskVal, stats); - - // 3. Use the sample median as the first estimator of the mean X. - clippedMean = stats->sampleMean; - - // 4. Use the sample stdev as the first estimator of the mean stdev. - clippedStdev = stats->sampleStdev; - - // Must save the old sampleMean and sampleStdev since the following code - // block overwrites them. - oldStanMean = stats->sampleMean; - oldStanStdev = stats->sampleStdev; - - // 5. Repeat N times: - for (i=0;iclipIter;i++) { - for (j=0;jn;j++) { - // a) Exclude all values x_i for which |x_i - x| > K * stdev - if (fabs(myVector->data.F32[j] - clippedMean) > - (stats->clipSigma * clippedStdev)) { - tmpMask->data.U8[i] = 0xff; - } - // b) compute new mean and stdev - p_psVectorSampleMedian(myVector, tmpMask, maskVal, stats); - p_psVectorSampleStdev(myVector, tmpMask, maskVal, stats); - - // c) Use the new mean for x - clippedMean = stats->sampleMean; - - // d) Use the new stdev for stdev - clippedStdev = stats->sampleStdev; - } - } - stats->sampleMean = oldStanMean; - stats->sampleStdev= oldStanStdev; - - // 7. The last calcuated value of x is the cliped mean. - if (stats->options & PS_STAT_CLIPPED_MEAN) { - stats->clippedMean = clippedMean; - } - - // 8. The last calcuated value of stdev is the cliped stdev. - if (stats->options & PS_STAT_CLIPPED_STDEV) { - stats->clippedStdev = clippedStdev; - } - - psVectorFree(tmpMask); }