Free cookie consent management tool by TermsFeed Policy Generator

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

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

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

File:
1 edited

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

 r13025 /// public static double Median(this IEnumerable values) { // iterate only once // See unit tests for comparison with naive implementation return Quantile(values, 0.5); } /// /// Calculates the alpha-quantile element of the enumeration. /// /// /// public static double Quantile(this IEnumerable 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; } } /// /// Calculates the alpha-quantile element of the enumeration. /// /// /// public static double Quantile(this IEnumerable 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. 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. } } } /// /// Calculates the pth percentile of the values. /// public static double Percentile(this IEnumerable values, double p) { // iterate only once
Note: See TracChangeset for help on using the changeset viewer.