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

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/sources/HeuristicLab.Common/3.3/EnumerableStatisticExtensions.cs
r13025 r13033 33 33 /// <returns></returns> 34 34 public static double Median(this IEnumerable<double> values) { 35 // iterate only once 35 // See unit tests for comparison with naive implementation 36 return Quantile(values, 0.5); 37 } 38 39 /// <summary> 40 /// Calculates the alphaquantile 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 36 46 double[] valuesArr = values.ToArray(); 37 47 int n = valuesArr.Length; 38 48 if (n == 0) throw new InvalidOperationException("Enumeration contains no elements."); 39 49 40 Array.Sort(valuesArr); 41 42 // return the middle element (if n is uneven) or the average of the two middle elements if n is even. 43 if (n % 2 == 1) { 44 return valuesArr[n / 2]; 45 } else { 46 return (valuesArr[(n / 2)  1] + valuesArr[n / 2]) / 2.0; 47 } 48 } 49 50 /// <summary> 51 /// Calculates the alphaquantile element of the enumeration. 52 /// </summary> 53 /// <param name="values"></param> 54 /// <returns></returns> 55 public static double Quantile(this IEnumerable<double> values, double alpha) { 56 Contract.Assert(alpha > 0 && alpha < 1); 57 // iterate only once 58 double[] valuesArr = values.ToArray(); 59 int n = valuesArr.Length; 60 if (n == 0) throw new InvalidOperationException("Enumeration contains no elements."); 61 62 Array.Sort(valuesArr); 63 // starts at 0 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." 64 53 65 54 // 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); 70 59 if (isInteger) { 71 return 0.5 * ( valuesArr[(int)pos  1] + valuesArr[(int)pos]);60 return 0.5 * (Select((int)pos  1, valuesArr) + Select((int)pos, valuesArr)); 72 61 } else { 73 return valuesArr[(int)Math.Ceiling(pos)  1]; 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..n1] returns an array value from array arr[0..n1] 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..k1] (in arbitrary order) and all 70 // larger elements in arr[k+1..n1] (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 } 74 141 } 75 142 } … … 133 200 /// <summary> 134 201 /// Calculates the pth percentile of the values. 135 /// </summary 202 /// </summary> 136 203 public static double Percentile(this IEnumerable<double> values, double p) { 137 204 // iterate only once
Note: See TracChangeset
for help on using the changeset viewer.