Changeset 13033 for trunk/sources/HeuristicLab.Common/3.3/EnumerableStatisticExtensions.cs

Timestamp:

10/19/15 16:13:14 (8 years ago)

Author:

gkronber

Message:

#2418: implemented O(n) algorithm for median (and general quantiles) determination.

File:

• trunk/sources/HeuristicLab.Common/3.3/EnumerableStatisticExtensions.cs (modified) (3 diffs)

trunk/sources/HeuristicLab.Common/3.3/EnumerableStatisticExtensions.cs

@@ -33 +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];
-      } else {
-        return (valuesArr[(n / 2) - 1] + valuesArr[n / 2]) / 2.0;
-      }
-    }
-
-    /// <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) {
-      Contract.Assert(alpha > 0 && alpha < 1);
-      // iterate only once 
-      double[] valuesArr = values.ToArray();
-      int n = valuesArr.Length;
-      if (n == 0) throw new InvalidOperationException("Enumeration contains no elements.");
-
-      Array.Sort(valuesArr);
-      //  starts at 0
+      // "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.
@@ -69 +58 @@
       bool isInteger = Math.Round(pos).IsAlmost(pos);
       if (isInteger) {
-        return 0.5 * (valuesArr[(int)pos - 1] + valuesArr[(int)pos]);
+        return 0.5 * (Select((int)pos - 1, valuesArr) + Select((int)pos, valuesArr));
       } else {
-        return valuesArr[(int)Math.Ceiling(pos) - 1];
+        return Select((int)Math.Ceiling(pos) - 1, valuesArr);
+      }
+    }
+
+    // Numerical Recipes in C++, §8.5 Selecting the Mth Largest, O(n)
+    // Giben k in [0..n-1] returns an array value from array arr[0..n-1] such that k array values are 
+    // lee than or equal to the one returned. The input array will be rearranged to hav 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).
+    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.
+        }
       }
     }
@@ -133 +200 @@
     /// <summary>
     /// Calculates the pth percentile of the values.
-    /// </summary
+    /// </summary>
     public static double Percentile(this IEnumerable<double> values, double p) {
       // iterate only once