Changeset 14029 for branches/crossvalidation-2434/HeuristicLab.Common/3.3/EnumerableStatisticExtensions.cs

Timestamp:

07/08/16 14:40:02 (8 years ago)

Author:

gkronber

Message:

#2434: merged trunk changes r12934:14026 from trunk to branch

Location:

branches/crossvalidation-2434

Files:

• . (modified) (1 prop)
• HeuristicLab.Common/3.3/EnumerableStatisticExtensions.cs (modified) (5 diffs)

branches/crossvalidation-2434/HeuristicLab.Common/3.3/EnumerableStatisticExtensions.cs

@@ -22 +22 @@
 using System;
 using System.Collections.Generic;
+using System.Diagnostics.Contracts;
 using System.Linq;
 
@@ -32 +33 @@
     /// <returns></returns>
     public static double Median(this IEnumerable<double> values) {
-      // iterate only once 
+      // See unit tests for comparison with naive implementation
+      return Quantile(values, 0.5);
+    }
+
+    /// <summary>
+    /// Calculates the alpha-quantile element of the enumeration.
+    /// </summary>
+    /// <param name="values"></param>
+    /// <returns></returns>
+    public static double Quantile(this IEnumerable<double> values, double alpha) {
+      // See unit tests for comparison with naive implementation
       double[] valuesArr = values.ToArray();
       int n = valuesArr.Length;
       if (n == 0) throw new InvalidOperationException("Enumeration contains no elements.");
 
-      Array.Sort(valuesArr);
-
-      // return the middle element (if n is uneven) or the average of the two middle elements if n is even.
-      if (n % 2 == 1) {
-        return valuesArr[n / 2];
+      // "When N is even, statistics books define the median as the arithmetic mean of the elements k = N/2 
+      // and k = N/2 + 1 (that is, N/2 from the bottom and N/2 from the top). 
+      // If you accept such pedantry, you must perform two separate selections to find these elements."
+
+      // return the element at Math.Ceiling (if n*alpha is fractional) or the average of two elements if n*alpha is integer.
+      var pos = n * alpha;
+      Contract.Assert(pos >= 0);
+      Contract.Assert(pos < n);
+      bool isInteger = Math.Round(pos).IsAlmost(pos);
+      if (isInteger) {
+        return 0.5 * (Select((int)pos - 1, valuesArr) + Select((int)pos, valuesArr));
       } else {
-        return (valuesArr[(n / 2) - 1] + valuesArr[n / 2]) / 2.0;
+        return Select((int)Math.Ceiling(pos) - 1, valuesArr);
+      }
+    }
+
+    // Numerical Recipes in C++, §8.5 Selecting the Mth Largest, O(n)
+    // Given k in [0..n-1] returns an array value from array arr[0..n-1] such that k array values are 
+    // less than or equal to the one returned. The input array will be rearranged to have this value in 
+    // location arr[k], with all smaller elements moved to arr[0..k-1] (in arbitrary order) and all 
+    // larger elements in arr[k+1..n-1] (also in arbitrary order).
+    // 
+    // Could be changed to Select<T> where T is IComparable but in this case is significantly slower for double values
+    private static double Select(int k, double[] arr) {
+      Contract.Assert(arr.GetLowerBound(0) == 0);
+      Contract.Assert(k >= 0 && k < arr.Length);
+      int i, ir, j, l, mid, n = arr.Length;
+      double a;
+      l = 0;
+      ir = n - 1;
+      for (; ; ) {
+        if (ir <= l + 1) {
+          // Active partition contains 1 or 2 elements.
+          if (ir == l + 1 && arr[ir] < arr[l]) {
+            // if (ir == l + 1 && arr[ir].CompareTo(arr[l]) < 0) {
+            // Case of 2 elements.
+            // SWAP(arr[l], arr[ir]);
+            double temp = arr[l];
+            arr[l] = arr[ir];
+            arr[ir] = temp;
+          }
+          return arr[k];
+        } else {
+          mid = (l + ir) >> 1; // Choose median of left, center, and right elements
+          {
+            // SWAP(arr[mid], arr[l + 1]); // as partitioning element a. Also
+            double temp = arr[mid];
+            arr[mid] = arr[l + 1];
+            arr[l + 1] = temp;
+          }
+
+          if (arr[l] > arr[ir]) {
+            // if (arr[l].CompareTo(arr[ir]) > 0) {  // rearrange so that arr[l] arr[ir] <= arr[l+1],
+            // SWAP(arr[l], arr[ir]); . arr[ir] >= arr[l+1]
+            double temp = arr[l];
+            arr[l] = arr[ir];
+            arr[ir] = temp;
+          }
+
+          if (arr[l + 1] > arr[ir]) {
+            // if (arr[l + 1].CompareTo(arr[ir]) > 0) {
+            // SWAP(arr[l + 1], arr[ir]);
+            double temp = arr[l + 1];
+            arr[l + 1] = arr[ir];
+            arr[ir] = temp;
+          }
+          if (arr[l] > arr[l + 1]) {
+            //if (arr[l].CompareTo(arr[l + 1]) > 0) {
+            // SWAP(arr[l], arr[l + 1]);
+            double temp = arr[l];
+            arr[l] = arr[l + 1];
+            arr[l + 1] = temp;
+
+          }
+          i = l + 1; // Initialize pointers for partitioning.
+          j = ir;
+          a = arr[l + 1]; // Partitioning element.
+          for (; ; ) { // Beginning of innermost loop.
+            do i++; while (arr[i] < a /* arr[i].CompareTo(a) < 0 */); // Scan up to find element > a.
+            do j--; while (arr[j] > a /* arr[j].CompareTo(a) > 0 */); // Scan down to find element < a.
+            if (j < i) break; // Pointers crossed. Partitioning complete.
+            {
+              // SWAP(arr[i], arr[j]);
+              double temp = arr[i];
+              arr[i] = arr[j];
+              arr[j] = temp;
+            }
+          } // End of innermost loop.
+          arr[l + 1] = arr[j]; // Insert partitioning element.
+          arr[j] = a;
+          if (j >= k) ir = j - 1; // Keep active the partition that contains the
+          if (j <= k) l = i; // kth element.
+        }
       }
     }
@@ -67 +164 @@
 
     /// <summary>
-    /// Calculates the standard deviation of values.
+    /// Calculates the sample standard deviation of values.
     /// </summary>
     /// <param name="values"></param>
@@ -76 +173 @@
 
     /// <summary>
-    /// Calculates the variance of values. (sum (x - x_mean)² / n)
+    /// Calculates the population standard deviation of values.
+    /// </summary>
+    /// <param name="values"></param>
+    /// <returns></returns>
+    public static double StandardDeviationPop(this IEnumerable<double> values) {
+      return Math.Sqrt(VariancePop(values));
+    }
+
+    /// <summary>
+    /// Calculates the sample variance of values. (sum (x - x_mean)² / (n-1))
     /// </summary>
     /// <param name="values"></param>
     /// <returns></returns>
     public static double Variance(this IEnumerable<double> values) {
+      return Variance(values, true);
+    }
+
+    /// <summary>
+    /// Calculates the population variance of values. (sum (x - x_mean)² / n)
+    /// </summary>
+    /// <param name="values"></param>
+    /// <returns></returns>
+    public static double VariancePop(this IEnumerable<double> values) {
+      return Variance(values, false);
+    }
+
+    private static double Variance(IEnumerable<double> values, bool sampleVariance) {
       int m_n = 0;
       double m_oldM = 0.0;
@@ -100 +219 @@
         }
       }
-      return ((m_n > 1) ? m_newS / (m_n - 1) : 0.0);
-    }
-
-    /// <summary>
-    /// Calculates the pth percentile of the values.
-    /// </summary
-    public static double Percentile(this IEnumerable<double> values, double p) {
-      // iterate only once 
-      double[] valuesArr = values.ToArray();
-      int n = valuesArr.Length;
-      if (n == 0) throw new InvalidOperationException("Enumeration contains no elements.");
-      if (n == 1) return values.ElementAt(0);
-
-      if (p.IsAlmost(0.0)) return valuesArr[0];
-      if (p.IsAlmost(1.0)) return valuesArr[n - 1];
-
-      double t = p * (n - 1);
-      int index = (int)Math.Floor(t);
-      double percentage = t - index;
-      return valuesArr[index] * (1 - percentage) + valuesArr[index + 1] * percentage;
+
+      if (m_n == 0) return double.NaN;
+      if (m_n == 1) return 0.0;
+
+      if (sampleVariance) return m_newS / (m_n - 1);
+      else return m_newS / m_n;
     }