Free cookie consent management tool by TermsFeed Policy Generator

Ignore:
Timestamp:
10/19/15 16:13:14 (9 years ago)
Author:
gkronber
Message:

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

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/sources/HeuristicLab.Tests/HeuristicLab.Common-3.3/EnumerableStatisticExtensions.cs

    r13025 r13033  
    2020#endregion
    2121
     22using System;
     23using System.Collections.Generic;
     24using System.Diagnostics;
     25using System.Diagnostics.Contracts;
    2226using System.Linq;
    2327using HeuristicLab.Common;
     28using HeuristicLab.Random;
    2429using Microsoft.VisualStudio.TestTools.UnitTesting;
    2530
     
    3136    [TestProperty("Time", "short")]
    3237    public void QuantileTest() {
     38      {
     39        Assert.AreEqual(2.0, new double[] { 2.0 }.Quantile(0.5));
     40        Assert.AreEqual(2.0, new double[] { 2.0 }.Quantile(0.01));
     41        Assert.AreEqual(2.0, new double[] { 2.0 }.Quantile(0.99));
     42      }
     43
     44      {
     45        Assert.AreEqual(1.5, new double[] { 1.0, 2.0 }.Median());
     46        Assert.AreEqual(2.0, new double[] { 1.0, 2.0 }.Quantile(0.99));
     47        Assert.AreEqual(1.0, new double[] { 1.0, 2.0 }.Quantile(0.01));
     48      }
     49      {
     50        Assert.AreEqual(2.0, new double[] { 3.0, 1.0, 2.0 }.Median());
     51        Assert.AreEqual(3.0, new double[] { 3.0, 1.0, 2.0 }.Quantile(0.99));
     52        Assert.AreEqual(1.0, new double[] { 3.0, 1.0, 2.0 }.Quantile(0.01));
     53      }
     54
     55
    3356      var xs = new double[] { 1, 1, 1, 3, 4, 7, 9, 11, 13, 13 };
    3457      {
    35         var q = xs.Quantile(0.3);
    36         Assert.AreEqual(q, 2.0, 1E-6);
     58        var q0 = Quantile(xs, 0.3); // naive implementation using sorting
     59        Assert.AreEqual(q0, 2.0, 1E-6);
     60
     61        var q1 = xs.Quantile(0.3); // using select
     62        Assert.AreEqual(q1, 2.0, 1E-6);
    3763      }
    3864      {
    39         var q = xs.Quantile(0.75);
    40         Assert.AreEqual(q, 11.0, 1E-6);
     65        var q0 = Quantile(xs, 0.75); // naive implementation using sorting
     66        Assert.AreEqual(q0, 11.0, 1E-6);
     67
     68        var q1 = xs.Quantile(0.75); // using select
     69        Assert.AreEqual(q1, 11.0, 1E-6);
    4170      }
    4271      // quantile = 0.5 is equivalent to median
    4372      {
    4473        // even number of elements
    45         Assert.AreEqual(xs.Quantile(0.5), xs.Median(), 1E-6);
     74        var expected = Median(xs);
     75        Assert.AreEqual(expected, Quantile(xs, 0.5), 1E-6); // using sorting
     76        Assert.AreEqual(expected, xs.Quantile(0.5), 1E-6); // using select
    4677      }
    4778      {
    4879        // odd number of elements
    49         Assert.AreEqual(xs.Take(9).Quantile(0.5), xs.Take(9).Median(), 1E-6);
     80        var expected = Median(xs.Take(9));
     81        Assert.AreEqual(expected, Quantile(xs.Take(9), 0.5), 1E-6); // using sorting
     82        Assert.AreEqual(expected, xs.Take(9).Quantile(0.5), 1E-6); // using select
     83      }
     84
     85      // edge cases
     86      {
     87        try {
     88          new double[] { }.Quantile(0.5); // empty
     89          Assert.Fail("expected exception");
     90        }
     91        catch (Exception) {
     92        }
     93      }
     94      {
     95        try {
     96          Enumerable.Repeat(0.0, 10).Quantile(1.0); // alpha < 1
     97          Assert.Fail("expected exception");
     98        }
     99        catch (Exception) {
     100        }
     101      }
     102    }
     103
     104    [TestMethod]
     105    [TestCategory("General")]
     106    [TestProperty("Time", "medium")]
     107    public void QuantilePerformanceTest() {
     108      int n = 10;
     109      var sw0 = new Stopwatch();
     110      var sw1 = new Stopwatch();
     111      const int reps = 1000;
     112      while (n <= 1000000) {
     113        for (int i = 0; i < reps; i++) {
     114          var xs = RandomEnumerable.SampleRandomNumbers(0, 10000, n + 1).Select(x => (double)x).ToArray();
     115          sw0.Start();
     116          var q0 = Median(xs); // sorting
     117          sw0.Stop();
     118
     119
     120          sw1.Start();
     121          var q1 = xs.Median(); // selection
     122          sw1.Stop();
     123          Assert.AreEqual(q0, q1, 1E-9);
     124        }
     125        Console.WriteLine("{0,-10} {1,-10} {2,-10}", n, sw0.ElapsedMilliseconds, sw1.ElapsedMilliseconds);
     126
     127        n = n * 10;
     128      }
     129    }
     130
     131
     132    // straight forward implementation of median function (using sorting)
     133    private static double Median(IEnumerable<double> values) {
     134      // iterate only once
     135      double[] valuesArr = values.ToArray();
     136      int n = valuesArr.Length;
     137      if (n == 0) throw new InvalidOperationException("Enumeration contains no elements.");
     138
     139      Array.Sort(valuesArr);
     140
     141      // return the middle element (if n is uneven) or the average of the two middle elements if n is even.
     142      if (n % 2 == 1) {
     143        return valuesArr[n / 2];
     144      } else {
     145        return (valuesArr[(n / 2) - 1] + valuesArr[n / 2]) / 2.0;
     146      }
     147    }
     148
     149    // straight forward implementation of quantile function (using sorting)
     150    private static double Quantile(IEnumerable<double> values, double alpha) {
     151      Contract.Assert(alpha > 0 && alpha < 1);
     152      // iterate only once
     153      double[] valuesArr = values.ToArray();
     154      int n = valuesArr.Length;
     155      if (n == 0) throw new InvalidOperationException("Enumeration contains no elements.");
     156
     157      Array.Sort(valuesArr);
     158      //  starts at 0
     159
     160      // return the element at Math.Ceiling (if n*alpha is fractional) or the average of two elements if n*alpha is integer.
     161      var pos = n * alpha;
     162      Contract.Assert(pos >= 0);
     163      Contract.Assert(pos < n);
     164      bool isInteger = Math.Round(pos).IsAlmost(pos);
     165      if (isInteger) {
     166        return 0.5 * (valuesArr[(int)pos - 1] + valuesArr[(int)pos]);
     167      } else {
     168        return valuesArr[(int)Math.Ceiling(pos) - 1];
    50169      }
    51170    }
Note: See TracChangeset for help on using the changeset viewer.