Changeset 13033

Timestamp:

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

Author:

gkronber

Message:

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

Location:

trunk/sources

Files:

• HeuristicLab.Common/3.3/EnumerableStatisticExtensions.cs (modified) (3 diffs)
• HeuristicLab.Tests/HeuristicLab.Common-3.3/EnumerableStatisticExtensions.cs (modified) (2 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 

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

@@ -20 +20 @@
 #endregion
 
+using System;
+using System.Collections.Generic;
+using System.Diagnostics;
+using System.Diagnostics.Contracts;
 using System.Linq;
 using HeuristicLab.Common;
+using HeuristicLab.Random;
 using Microsoft.VisualStudio.TestTools.UnitTesting;
 
@@ -31 +36 @@
     [TestProperty("Time", "short")]
     public void QuantileTest() {
+      {
+        Assert.AreEqual(2.0, new double[] { 2.0 }.Quantile(0.5));
+        Assert.AreEqual(2.0, new double[] { 2.0 }.Quantile(0.01));
+        Assert.AreEqual(2.0, new double[] { 2.0 }.Quantile(0.99));
+      }
+
+      {
+        Assert.AreEqual(1.5, new double[] { 1.0, 2.0 }.Median());
+        Assert.AreEqual(2.0, new double[] { 1.0, 2.0 }.Quantile(0.99));
+        Assert.AreEqual(1.0, new double[] { 1.0, 2.0 }.Quantile(0.01));
+      }
+      {
+        Assert.AreEqual(2.0, new double[] { 3.0, 1.0, 2.0 }.Median());
+        Assert.AreEqual(3.0, new double[] { 3.0, 1.0, 2.0 }.Quantile(0.99));
+        Assert.AreEqual(1.0, new double[] { 3.0, 1.0, 2.0 }.Quantile(0.01));
+      }
+
+
       var xs = new double[] { 1, 1, 1, 3, 4, 7, 9, 11, 13, 13 };
       {
-        var q = xs.Quantile(0.3);
-        Assert.AreEqual(q, 2.0, 1E-6);
+        var q0 = Quantile(xs, 0.3); // naive implementation using sorting
+        Assert.AreEqual(q0, 2.0, 1E-6);
+
+        var q1 = xs.Quantile(0.3); // using select
+        Assert.AreEqual(q1, 2.0, 1E-6);
       }
       {
-        var q = xs.Quantile(0.75);
-        Assert.AreEqual(q, 11.0, 1E-6);
+        var q0 = Quantile(xs, 0.75); // naive implementation using sorting
+        Assert.AreEqual(q0, 11.0, 1E-6);
+
+        var q1 = xs.Quantile(0.75); // using select
+        Assert.AreEqual(q1, 11.0, 1E-6);
       }
       // quantile = 0.5 is equivalent to median
       {
         // even number of elements
-        Assert.AreEqual(xs.Quantile(0.5), xs.Median(), 1E-6);
+        var expected = Median(xs);
+        Assert.AreEqual(expected, Quantile(xs, 0.5), 1E-6); // using sorting
+        Assert.AreEqual(expected, xs.Quantile(0.5), 1E-6); // using select
       }
       {
         // odd number of elements
-        Assert.AreEqual(xs.Take(9).Quantile(0.5), xs.Take(9).Median(), 1E-6);
+        var expected = Median(xs.Take(9));
+        Assert.AreEqual(expected, Quantile(xs.Take(9), 0.5), 1E-6); // using sorting
+        Assert.AreEqual(expected, xs.Take(9).Quantile(0.5), 1E-6); // using select
+      }
+
+      // edge cases
+      {
+        try {
+          new double[] { }.Quantile(0.5); // empty
+          Assert.Fail("expected exception");
+        }
+        catch (Exception) {
+        }
+      }
+      {
+        try {
+          Enumerable.Repeat(0.0, 10).Quantile(1.0); // alpha < 1 
+          Assert.Fail("expected exception");
+        }
+        catch (Exception) {
+        }
+      }
+    }
+
+    [TestMethod]
+    [TestCategory("General")]
+    [TestProperty("Time", "medium")]
+    public void QuantilePerformanceTest() {
+      int n = 10;
+      var sw0 = new Stopwatch();
+      var sw1 = new Stopwatch();
+      const int reps = 1000;
+      while (n <= 1000000) {
+        for (int i = 0; i < reps; i++) {
+          var xs = RandomEnumerable.SampleRandomNumbers(0, 10000, n + 1).Select(x => (double)x).ToArray();
+          sw0.Start();
+          var q0 = Median(xs); // sorting
+          sw0.Stop();
+
+
+          sw1.Start();
+          var q1 = xs.Median(); // selection
+          sw1.Stop();
+          Assert.AreEqual(q0, q1, 1E-9);
+        }
+        Console.WriteLine("{0,-10} {1,-10} {2,-10}", n, sw0.ElapsedMilliseconds, sw1.ElapsedMilliseconds);
+
+        n = n * 10;
+      }
+    }
+
+
+    // straight forward implementation of median function (using sorting)
+    private static double Median(IEnumerable<double> values) {
+      // iterate only once 
+      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;
+      }
+    }
+
+    // straight forward implementation of quantile function (using sorting)
+    private static double Quantile(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
+
+      // 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 * (valuesArr[(int)pos - 1] + valuesArr[(int)pos]);
+      } else {
+        return valuesArr[(int)Math.Ceiling(pos) - 1];
       }
     }