source: trunk/sources/HeuristicLab.Common/3.3/EnumerableStatisticExtensions.cs @ 13034

Last change on this file since 13034 was 13034, checked in by gkronber, 4 years ago

#2491: removed separate (incorrect) implementation of percentile function

File size: 8.2 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2015 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
4 *
5 * This file is part of HeuristicLab.
6 *
7 * HeuristicLab is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * HeuristicLab is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
19 */
20#endregion
21
22using System;
23using System.Collections.Generic;
24using System.Diagnostics.Contracts;
25using System.Linq;
26
27namespace HeuristicLab.Common {
28  public static class EnumerableStatisticExtensions {
29    /// <summary>
30    /// Calculates the median element of the enumeration.
31    /// </summary>
32    /// <param name="values"></param>
33    /// <returns></returns>
34    public static double Median(this IEnumerable<double> values) {
35      // See unit tests for comparison with naive implementation
36      return Quantile(values, 0.5);
37    }
38
39    /// <summary>
40    /// Calculates the alpha-quantile element of the enumeration.
41    /// </summary>
42    /// <param name="values"></param>
43    /// <returns></returns>
44    public static double Quantile(this IEnumerable<double> values, double alpha) {
45      // See unit tests for comparison with naive implementation
46      double[] valuesArr = values.ToArray();
47      int n = valuesArr.Length;
48      if (n == 0) throw new InvalidOperationException("Enumeration contains no elements.");
49
50      // "When N is even, statistics books define the median as the arithmetic mean of the elements k = N/2
51      // and k = N/2 + 1 (that is, N/2 from the bottom and N/2 from the top).
52      // If you accept such pedantry, you must perform two separate selections to find these elements."
53
54      // return the element at Math.Ceiling (if n*alpha is fractional) or the average of two elements if n*alpha is integer.
55      var pos = n * alpha;
56      Contract.Assert(pos >= 0);
57      Contract.Assert(pos < n);
58      bool isInteger = Math.Round(pos).IsAlmost(pos);
59      if (isInteger) {
60        return 0.5 * (Select((int)pos - 1, valuesArr) + Select((int)pos, valuesArr));
61      } else {
62        return Select((int)Math.Ceiling(pos) - 1, valuesArr);
63      }
64    }
65
66    // Numerical Recipes in C++, §8.5 Selecting the Mth Largest, O(n)
67    // Giben k in [0..n-1] returns an array value from array arr[0..n-1] such that k array values are
68    // lee than or equal to the one returned. The input array will be rearranged to hav this value in
69    // location arr[k], with all smaller elements moved to arr[0..k-1] (in arbitrary order) and all
70    // larger elements in arr[k+1..n-1] (also in arbitrary order).
71    private static double Select(int k, double[] arr) {
72      Contract.Assert(arr.GetLowerBound(0) == 0);
73      Contract.Assert(k >= 0 && k < arr.Length);
74      int i, ir, j, l, mid, n = arr.Length;
75      double a;
76      l = 0;
77      ir = n - 1;
78      for (; ; ) {
79        if (ir <= l + 1) {
80          // Active partition contains 1 or 2 elements.
81          if (ir == l + 1 && arr[ir] < arr[l]) {
82            // if (ir == l + 1 && arr[ir].CompareTo(arr[l]) < 0) {
83            // Case of 2 elements.
84            // SWAP(arr[l], arr[ir]);
85            double temp = arr[l];
86            arr[l] = arr[ir];
87            arr[ir] = temp;
88          }
89          return arr[k];
90        } else {
91          mid = (l + ir) >> 1; // Choose median of left, center, and right elements
92          {
93            // SWAP(arr[mid], arr[l + 1]); // as partitioning element a. Also
94            double temp = arr[mid];
95            arr[mid] = arr[l + 1];
96            arr[l + 1] = temp;
97          }
98
99          if (arr[l] > arr[ir]) {
100            // if (arr[l].CompareTo(arr[ir]) > 0) {  // rearrange so that arr[l] arr[ir] <= arr[l+1],
101            // SWAP(arr[l], arr[ir]); . arr[ir] >= arr[l+1]
102            double temp = arr[l];
103            arr[l] = arr[ir];
104            arr[ir] = temp;
105          }
106
107          if (arr[l + 1] > arr[ir]) {
108            // if (arr[l + 1].CompareTo(arr[ir]) > 0) {
109            // SWAP(arr[l + 1], arr[ir]);
110            double temp = arr[l + 1];
111            arr[l + 1] = arr[ir];
112            arr[ir] = temp;
113          }
114          if (arr[l] > arr[l + 1]) {
115          //if (arr[l].CompareTo(arr[l + 1]) > 0) {
116            // SWAP(arr[l], arr[l + 1]);
117            double temp = arr[l];
118            arr[l] = arr[l + 1];
119            arr[l + 1] = temp;
120
121          }
122          i = l + 1; // Initialize pointers for partitioning.
123          j = ir;
124          a = arr[l + 1]; // Partitioning element.
125          for (; ; ) { // Beginning of innermost loop.
126            do i++; while (arr[i] < a /* arr[i].CompareTo(a) < 0 */); // Scan up to find element > a.
127            do j--; while (arr[j] > a /* arr[j].CompareTo(a) > 0 */); // Scan down to find element < a.
128            if (j < i) break; // Pointers crossed. Partitioning complete.
129            {
130              // SWAP(arr[i], arr[j]);
131              double temp = arr[i];
132              arr[i] = arr[j];
133              arr[j] = temp;
134            }
135          } // End of innermost loop.
136          arr[l + 1] = arr[j]; // Insert partitioning element.
137          arr[j] = a;
138          if (j >= k) ir = j - 1; // Keep active the partition that contains the
139          if (j <= k) l = i; // kth element.
140        }
141      }
142    }
143
144    /// <summary>
145    /// Calculates the range (max - min) of the enumeration.
146    /// </summary>
147    /// <param name="values"></param>
148    /// <returns></returns>
149    public static double Range(this IEnumerable<double> values) {
150      double min = double.PositiveInfinity;
151      double max = double.NegativeInfinity;
152      int i = 0;
153      foreach (var e in values) {
154        if (min > e) min = e;
155        if (max < e) max = e;
156        i++;
157      }
158      if (i < 1) throw new ArgumentException("The enumerable must contain at least two elements", "values");
159      return max - min;
160    }
161
162
163    /// <summary>
164    /// Calculates the standard deviation of values.
165    /// </summary>
166    /// <param name="values"></param>
167    /// <returns></returns>
168    public static double StandardDeviation(this IEnumerable<double> values) {
169      return Math.Sqrt(Variance(values));
170    }
171
172    /// <summary>
173    /// Calculates the variance of values. (sum (x - x_mean)² / n)
174    /// </summary>
175    /// <param name="values"></param>
176    /// <returns></returns>
177    public static double Variance(this IEnumerable<double> values) {
178      int m_n = 0;
179      double m_oldM = 0.0;
180      double m_newM = 0.0;
181      double m_oldS = 0.0;
182      double m_newS = 0.0;
183      foreach (double x in values) {
184        m_n++;
185        if (m_n == 1) {
186          m_oldM = m_newM = x;
187          m_oldS = 0.0;
188        } else {
189          m_newM = m_oldM + (x - m_oldM) / m_n;
190          m_newS = m_oldS + (x - m_oldM) * (x - m_newM);
191
192          // set up for next iteration
193          m_oldM = m_newM;
194          m_oldS = m_newS;
195        }
196      }
197      return ((m_n > 1) ? m_newS / (m_n - 1) : 0.0);
198    }
199
200    /// <summary>
201    /// Calculates the pth percentile of the values.
202    /// </summary>
203    public static double Percentile(this IEnumerable<double> values, double p) {
204      // percentiles are actually quantiles where alpha is constrained to integer percentage values from 1% to 99%
205      return Quantile(values, p);
206    }
207
208    public static IEnumerable<double> LimitToRange(this IEnumerable<double> values, double min, double max) {
209      if (min > max) throw new ArgumentException(string.Format("Minimum {0} is larger than maximum {1}.", min, max));
210      foreach (var x in values) {
211        if (double.IsNaN(x)) yield return (max + min) / 2.0;
212        else if (x < min) yield return min;
213        else if (x > max) yield return max;
214        else yield return x;
215      }
216    }
217  }
218}
Note: See TracBrowser for help on using the repository browser.