Changeset 14428

Timestamp:

Aug 8, 2007, 9:34:58 AM (19 years ago)

Author:

Paul Price

Message:

Adding higher-order moments (skewness and kurtosis).

File:

• trunk/psLib/src/math/psStats.c (modified) (10 diffs)

trunk/psLib/src/math/psStats.c

@@ -13 +13 @@
  * 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 +382 @@
 
 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 +461 @@
     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 +1753 @@
     stats->sampleUQ = NAN;
     stats->sampleLQ = NAN;
+    stats->sampleSkewness = NAN;
+    stats->sampleKurtosis = NAN;
     stats->robustMedian = NAN;
     stats->robustStdev = NAN;
@@ -1808 +1880 @@
     }
 
+    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 +1957 @@
     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 +1965 @@
     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 +2002 @@
     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 +2056 @@
       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 +2089 @@
       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;