Index: trunk/psLib/src/math/psStats.c =================================================================== --- trunk/psLib/src/math/psStats.c (revision 13991) +++ trunk/psLib/src/math/psStats.c (revision 14428) @@ -13,6 +13,6 @@ * use ->min and ->max (PS_STAT_USE_RANGE) * - * @version $Revision: 1.213 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-06-30 00:35:38 $ + * @version $Revision: 1.214 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-08-08 19:34:58 $ * * Copyright 2006 IfA, University of Hawaii @@ -382,18 +382,12 @@ 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 - +much easier to debug. *****************************************************************************/ - static bool vectorSampleStdev(const psVector* myVector, - const psVector* errors, - const psVector* maskVector, - psMaskType maskVal, - psStats* stats) + +static bool vectorSampleStdev(const psVector* myVector, + const psVector* errors, + const psVector* maskVector, + psMaskType maskVal, + psStats* stats) { psTrace(TRACE, 4, "---- %s() begin ----\n", __func__); @@ -467,4 +461,80 @@ psTrace(TRACE, 4, "---- %s() end ----\n", __func__); + return true; +} + +static bool vectorSampleMoments(const psVector* myVector, + const psVector* maskVector, + psMaskType maskVal, + psStats* stats) +{ + psTrace(TRACE, 4, "---- %s() begin ----\n", __func__); + + // This procedure requires the mean and standard deviation + if (!(stats->results & PS_STAT_SAMPLE_MEAN)) { + vectorSampleMean(myVector, NULL, maskVector, maskVal, stats); + } + if (isnan(stats->sampleMean)) { + psTrace(TRACE, 5, "WARNING: vectorSampleMoments(): sample mean is NAN.\n"); + goto SAMPLE_MOMENTS_BAD; + } + if (!(stats->results & PS_STAT_SAMPLE_STDEV)) { + vectorSampleMean(myVector, NULL, maskVector, maskVal, stats); + } + if (isnan(stats->sampleStdev) || stats->sampleStdev == 0.0) { + psTrace(TRACE, 5, "WARNING: vectorSampleMoments(): sample mean is NAN or 0.\n"); + goto SAMPLE_MOMENTS_BAD; + } + + psF32 *data = myVector->data.F32; // Dereference + psU8 *maskData = (maskVector == NULL) ? NULL : maskVector->data.U8; + bool useRange = stats->options & PS_STAT_USE_RANGE; + + // Accumulate the sums + double mean = stats->sampleMean; // The mean + double sum3 = 0.0; // Sum of the cubes of the differences + double sum4 = 0.0; // Sum of the fourth powers of the differences + long count = 0; // Number of data points being used + 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; // Difference from the mean + double temp; // Temporary variable for accumulating + + sum3 += temp = PS_SQR(diff); + sum4 += temp * diff; + + count++; + } + + assert(count > 1); // It should be, because we have a mean and standard deviation + + double stdev = stats->sampleStdev; // Standard deviation + double variance = PS_SQR(stdev); // Variance + + // Formula for skewness and kurtosis from Numerical Recipes in C, p 613. + // Note that we are defining the kurtosis relative to a normal distribution (hence the -3.0). + stats->sampleSkewness = sum3 / (count * variance * stdev); + stats->sampleKurtosis = sum4 / (count * PS_SQR(variance)) - 3.0; + stats->results |= PS_STAT_SAMPLE_SKEWNESS | PS_STAT_SAMPLE_KURTOSIS; + + psTrace(TRACE, 4, "---- %s() end ----\n", __func__); + + return true; + + SAMPLE_MOMENTS_BAD: + // stats->sampleStdev has already been set + stats->sampleSkewness = NAN; + stats->sampleKurtosis = NAN; + stats->results |= PS_STAT_SAMPLE_SKEWNESS | PS_STAT_SAMPLE_KURTOSIS; return true; } @@ -1683,4 +1753,6 @@ stats->sampleUQ = NAN; stats->sampleLQ = NAN; + stats->sampleSkewness = NAN; + stats->sampleKurtosis = NAN; stats->robustMedian = NAN; stats->robustStdev = NAN; @@ -1808,4 +1880,11 @@ } + if (stats->options & (PS_STAT_SAMPLE_SKEWNESS | PS_STAT_SAMPLE_KURTOSIS)) { + if (!vectorSampleMoments(inF32, maskU8, maskVal, stats)) { + psError(PS_ERR_UNKNOWN, false, "Failed to calculate sample moments"); + status &= false; + } + } + // ************************************************************************ if (stats->options & (PS_STAT_MAX | PS_STAT_MIN)) { @@ -1878,4 +1957,6 @@ READ_STAT("MEAN", PS_STAT_SAMPLE_MEAN); READ_STAT("STDEV", PS_STAT_SAMPLE_STDEV); + READ_STAT("SKEWNESS", PS_STAT_SAMPLE_SKEWNESS); + READ_STAT("KURTOSIS", PS_STAT_SAMPLE_KURTOSIS); READ_STAT("MEDIAN", PS_STAT_SAMPLE_MEDIAN); READ_STAT("QUARTILE", PS_STAT_SAMPLE_QUARTILE); @@ -1884,4 +1965,6 @@ READ_STAT("SAMPLE_MEDIAN", PS_STAT_SAMPLE_MEDIAN); READ_STAT("SAMPLE_QUARTILE", PS_STAT_SAMPLE_QUARTILE); + READ_STAT("SAMPLE_SKEWNESS", PS_STAT_SAMPLE_SKEWNESS); + READ_STAT("SAMPLE_KURTOSIS", PS_STAT_SAMPLE_KURTOSIS); READ_STAT("ROBUST", PS_STAT_ROBUST_MEDIAN); READ_STAT("ROBUST_MEDIAN", PS_STAT_ROBUST_MEDIAN); @@ -1919,4 +2002,6 @@ WRITE_STAT("SAMPLE_MEDIAN", PS_STAT_SAMPLE_MEDIAN); WRITE_STAT("SAMPLE_QUARTILE", PS_STAT_SAMPLE_QUARTILE); + WRITE_STAT("SAMPLE_SKEWNESS", PS_STAT_SAMPLE_SKEWNESS); + WRITE_STAT("SAMPLE_KURTOSIS", PS_STAT_SAMPLE_KURTOSIS); WRITE_STAT("ROBUST_MEDIAN", PS_STAT_ROBUST_MEDIAN); WRITE_STAT("ROBUST_STDEV", PS_STAT_ROBUST_STDEV); @@ -1971,4 +2056,6 @@ case PS_STAT_SAMPLE_STDEV: case PS_STAT_SAMPLE_QUARTILE: + case PS_STAT_SAMPLE_SKEWNESS: + case PS_STAT_SAMPLE_KURTOSIS: case PS_STAT_ROBUST_MEDIAN: case PS_STAT_ROBUST_STDEV: @@ -2002,4 +2089,8 @@ case PS_STAT_SAMPLE_STDEV: return stats->sampleStdev; + case PS_STAT_SAMPLE_SKEWNESS: + return stats->sampleSkewness; + case PS_STAT_SAMPLE_KURTOSIS: + return stats->sampleKurtosis; case PS_STAT_ROBUST_MEDIAN: return stats->robustMedian;