From 0b5f954739944c523a099e494ace612517490948 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sun, 29 Dec 2013 16:39:25 +0100 Subject: [PATCH] Statistics: Ranks --- src/Numerics/Numerics.csproj | 1 + src/Numerics/Sorting.cs | 265 ++++++++++-------- src/Numerics/Statistics/ArrayStatistics.cs | 210 +++++++++----- src/Numerics/Statistics/Correlation.cs | 37 +-- src/Numerics/Statistics/RankDefinition.cs | 49 ++++ .../Statistics/SortedArrayStatistics.cs | 196 +++++++++---- src/Numerics/Statistics/Statistics.cs | 26 ++ src/UnitTests/ExcelTests.cs | 7 + .../StatisticsTests/StatisticsTests.cs | 121 ++++++++ 9 files changed, 638 insertions(+), 274 deletions(-) create mode 100644 src/Numerics/Statistics/RankDefinition.cs diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 624bd6de..31785ea8 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -202,6 +202,7 @@ + diff --git a/src/Numerics/Sorting.cs b/src/Numerics/Sorting.cs index e8e79f53..b9a813b2 100644 --- a/src/Numerics/Sorting.cs +++ b/src/Numerics/Sorting.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2010 Math.NET +// Copyright (c) 2009-2013 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -28,52 +28,16 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using System; +using System.Collections.Generic; + namespace MathNet.Numerics { - using System; - using System.Collections.Generic; - /// /// Sorting algorithms for single, tuple and triple lists. /// public static class Sorting { - /// - /// Sort a list of keys, in place using the quick sort algorithm. - /// - /// The type of elements stored in the list. - /// List to sort. - public static void Sort(IList keys) - { - Sort(keys, Comparer.Default); - } - - /// - /// Sort a list of keys and items with respect to the keys, in place using the quick sort algorithm. - /// - /// The type of elements stored in the key list. - /// The type of elements stored in the item list. - /// List to sort. - /// List to permute the same way as the key list. - public static void Sort(IList keys, IList items) - { - Sort(keys, items, Comparer.Default); - } - - /// - /// Sort a list of keys, items1 and items2 with respect to the keys, in place using the quick sort algorithm. - /// - /// The type of elements stored in the key list. - /// The type of elements stored in the first item list. - /// The type of elements stored in the second item list. - /// List to sort. - /// First list to permute the same way as the key list. - /// Second list to permute the same way as the key list. - public static void Sort(IList keys, IList items1, IList items2) - { - Sort(keys, items1, items2, Comparer.Default); - } - /// /// Sort a range of a list of keys, in place using the quick sort algorithm. /// @@ -92,16 +56,11 @@ namespace MathNet.Numerics /// The type of elements in the key list. /// List to sort. /// Comparison, defining the sort order. - public static void Sort(IList keys, IComparer comparer) + public static void Sort(IList keys, IComparer comparer = null) { - if (null == keys) - { - throw new ArgumentNullException("keys"); - } - if (null == comparer) { - throw new ArgumentNullException("comparer"); + comparer = Comparer.Default; } // basic cases @@ -148,21 +107,11 @@ namespace MathNet.Numerics /// List to sort. /// List to permute the same way as the key list. /// Comparison, defining the sort order. - public static void Sort(IList keys, IList items, IComparer comparer) + public static void Sort(IList keys, IList items, IComparer comparer = null) { - if (null == keys) - { - throw new ArgumentNullException("keys"); - } - - if (null == items) - { - throw new ArgumentNullException("items"); - } - if (null == comparer) { - throw new ArgumentNullException("comparer"); + comparer = Comparer.Default; } #if !PORTABLE @@ -190,31 +139,40 @@ namespace MathNet.Numerics /// First list to permute the same way as the key list. /// Second list to permute the same way as the key list. /// Comparison, defining the sort order. - public static void Sort( - IList keys, IList items1, IList items2, IComparer comparer) + public static void Sort(IList keys, IList items1, IList items2, IComparer comparer = null) { - if (null == keys) + if (null == comparer) { - throw new ArgumentNullException("keys"); + comparer = Comparer.Default; } - if (null == items1) - { - throw new ArgumentNullException("items1"); - } + // local sort implementation + QuickSort(keys, items1, items2, comparer, 0, keys.Count - 1); + } - if (null == items2) + /// + /// Sort a list of keys and items with respect to the keys, in place using the quick sort algorithm. + /// + /// The type of elements in the primary list. + /// The type of elements in the secondary list. + /// List to sort. + /// List to to sort on duplicate primaty items, and permute the same way as the key list. + /// Comparison, defining the primary sort order. + /// Comparison, defining the secondary sort order. + public static void SortAll(IList primary, IList secondary, IComparer primaryComparer = null, IComparer secondaryComparer = null) + { + if (null == primaryComparer) { - throw new ArgumentNullException("items2"); + primaryComparer = Comparer.Default; } - if (null == comparer) + if (null == secondaryComparer) { - throw new ArgumentNullException("comparer"); + secondaryComparer = Comparer.Default; } // local sort implementation - QuickSort(keys, items1, items2, comparer, 0, keys.Count - 1); + QuickSortAll(primary, secondary, primaryComparer, secondaryComparer, 0, primary.Count - 1); } /// @@ -225,18 +183,8 @@ namespace MathNet.Numerics /// The zero-based starting index of the range to sort. /// The length of the range to sort. /// Comparison, defining the sort order. - public static void Sort(IList keys, int index, int count, IComparer comparer) + public static void Sort(IList keys, int index, int count, IComparer comparer = null) { - if (null == keys) - { - throw new ArgumentNullException("keys"); - } - - if (null == comparer) - { - throw new ArgumentNullException("comparer"); - } - if (index < 0 || index >= keys.Count) { throw new ArgumentOutOfRangeException("index"); @@ -247,6 +195,11 @@ namespace MathNet.Numerics throw new ArgumentOutOfRangeException("count"); } + if (null == comparer) + { + comparer = Comparer.Default; + } + // basic cases if (count <= 1) { @@ -291,11 +244,7 @@ namespace MathNet.Numerics /// The method with which to compare two elements of the quick sort. /// The left boundary of the quick sort. /// The right boundary of the quick sort. - private static void QuickSort( - IList keys, - IComparer comparer, - int left, - int right) + static void QuickSort(IList keys, IComparer comparer, int left, int right) { do { @@ -346,10 +295,9 @@ namespace MathNet.Numerics a++; b--; - } - while (a <= b); + } while (a <= b); - // In order to limit the recusion depth to log(n), we sort the + // In order to limit the recusion depth to log(n), we sort the // shorter partition recusively and the longer partition iteratively. if ((b - left) <= (right - a)) { @@ -369,8 +317,7 @@ namespace MathNet.Numerics right = b; } - } - while (left < right); + } while (left < right); } /// @@ -383,12 +330,7 @@ namespace MathNet.Numerics /// The method with which to compare two elements of the quick sort. /// The left boundary of the quick sort. /// The right boundary of the quick sort. - private static void QuickSort( - IList keys, - IList items, - IComparer comparer, - int left, - int right) + static void QuickSort(IList keys, IList items, IComparer comparer, int left, int right) { do { @@ -443,10 +385,9 @@ namespace MathNet.Numerics a++; b--; - } - while (a <= b); + } while (a <= b); - // In order to limit the recusion depth to log(n), we sort the + // In order to limit the recusion depth to log(n), we sort the // shorter partition recusively and the longer partition iteratively. if ((b - left) <= (right - a)) { @@ -466,8 +407,7 @@ namespace MathNet.Numerics right = b; } - } - while (left < right); + } while (left < right); } /// @@ -482,13 +422,10 @@ namespace MathNet.Numerics /// The method with which to compare two elements of the quick sort. /// The left boundary of the quick sort. /// The right boundary of the quick sort. - private static void QuickSort( - IList keys, - IList items1, - IList items2, + static void QuickSort( + IList keys, IList items1, IList items2, IComparer comparer, - int left, - int right) + int left, int right) { do { @@ -547,10 +484,9 @@ namespace MathNet.Numerics a++; b--; - } - while (a <= b); + } while (a <= b); - // In order to limit the recusion depth to log(n), we sort the + // In order to limit the recusion depth to log(n), we sort the // shorter partition recusively and the longer partition iteratively. if ((b - left) <= (right - a)) { @@ -570,8 +506,107 @@ namespace MathNet.Numerics right = b; } - } - while (left < right); + } while (left < right); + } + + /// + /// Recursive implementation for an in place quick sort on the primary and then by the secondary list while reordering one secondary list accordingly. + /// + /// The type of the primary list. + /// The type of the secondary list. + /// The list which is sorted using quick sort. + /// The list which is sorted secondarily (on primary duplicates) and automatically reordered accordingly. + /// The method with which to compare two elements of the primary list. + /// The method with which to compare two elements of the secondary list. + /// The left boundary of the quick sort. + /// The right boundary of the quick sort. + static void QuickSortAll( + IList primary, IList secondary, + IComparer primaryComparer, IComparer secondaryComparer, + int left, int right) + { + do + { + // Pivoting + int a = left; + int b = right; + int p = a + ((b - a) >> 1); // midpoint + + int ap = primaryComparer.Compare(primary[a], primary[p]); + if (ap > 0 || ap == 0 && secondaryComparer.Compare(secondary[a], secondary[p]) > 0) + { + Swap(primary, a, p); + Swap(secondary, a, p); + } + + int ab = primaryComparer.Compare(primary[a], primary[b]); + if (ab > 0 || ab == 0 && secondaryComparer.Compare(secondary[a], secondary[b]) > 0) + { + Swap(primary, a, b); + Swap(secondary, a, b); + } + + int pb = primaryComparer.Compare(primary[p], primary[b]); + if (pb > 0 || pb == 0 && secondaryComparer.Compare(secondary[p], secondary[b]) > 0) + { + Swap(primary, p, b); + Swap(secondary, p, b); + } + + T1 pivot1 = primary[p]; + T2 pivot2 = secondary[p]; + + // Hoare Partitioning + do + { + int ax; + while ((ax = primaryComparer.Compare(primary[a], pivot1)) < 0 || ax == 0 && secondaryComparer.Compare(secondary[a], pivot2) < 0) + { + a++; + } + + int xb; + while ((xb = primaryComparer.Compare(pivot1, primary[b])) < 0 || xb == 0 && secondaryComparer.Compare(pivot2, secondary[b]) < 0) + { + b--; + } + + if (a > b) + { + break; + } + + if (a < b) + { + Swap(primary, a, b); + Swap(secondary, a, b); + } + + a++; + b--; + } while (a <= b); + + // In order to limit the recusion depth to log(n), we sort the + // shorter partition recusively and the longer partition iteratively. + if ((b - left) <= (right - a)) + { + if (left < b) + { + QuickSortAll(primary, secondary, primaryComparer, secondaryComparer, left, b); + } + + left = a; + } + else + { + if (a < right) + { + QuickSortAll(primary, secondary, primaryComparer, secondaryComparer, a, right); + } + + right = b; + } + } while (left < right); } /// diff --git a/src/Numerics/Statistics/ArrayStatistics.cs b/src/Numerics/Statistics/ArrayStatistics.cs index 2877f0c3..991b9ce7 100644 --- a/src/Numerics/Statistics/ArrayStatistics.cs +++ b/src/Numerics/Statistics/ArrayStatistics.cs @@ -43,7 +43,7 @@ namespace MathNet.Numerics.Statistics public static class ArrayStatistics { // TODO: Benchmark various options to find out the best approach (-> branch prediction) - // TODO: consider leveraging MKL + // TODO: consider leveraging MKL /// /// Returns the smallest value from the unsorted data array. @@ -257,7 +257,7 @@ namespace MathNet.Numerics.Statistics var covariance = 0.0; for (int i = 0; i < population1.Length; i++) { - covariance += (population1[i] - mean1) * (population2[i] - mean2); + covariance += (population1[i] - mean1)*(population2[i] - mean2); } return covariance/population1.Length; } @@ -299,7 +299,7 @@ namespace MathNet.Numerics.Statistics /// Percentile selector, between 0 and 100 (inclusive). public static double PercentileInplace(double[] data, int p) { - return QuantileInplace(data, p / 100d); + return QuantileInplace(data, p/100d); } /// @@ -368,7 +368,7 @@ namespace MathNet.Numerics.Statistics if (tau < 0d || tau > 1d || data.Length == 0) return double.NaN; double h = (data.Length + 1d/3d)*tau + 1d/3d; - var hf = (int) h; + var hf = (int)h; if (hf <= 0 || tau == 0d) { @@ -401,7 +401,7 @@ namespace MathNet.Numerics.Statistics { if (tau < 0d || tau > 1d || data.Length == 0) return double.NaN; - var x = a + (data.Length + b) * tau - 1; + var x = a + (data.Length + b)*tau - 1; #if PORTABLE var ip = (int)x; #else @@ -411,12 +411,12 @@ namespace MathNet.Numerics.Statistics if (Math.Abs(fp) < 1e-9) { - return SelectInplace(data, (int) ip); + return SelectInplace(data, (int)ip); } - var lower = SelectInplace(data, (int) Math.Floor(x)); - var upper = SelectInplace(data, (int) Math.Ceiling(x)); - return lower + (upper - lower) * (c + d * fp); + var lower = SelectInplace(data, (int)Math.Floor(x)); + var upper = SelectInplace(data, (int)Math.Ceiling(x)); + return lower + (upper - lower)*(c + d*fp); } /// @@ -438,68 +438,68 @@ namespace MathNet.Numerics.Statistics switch (definition) { case QuantileDefinition.R1: - { - double h = data.Length * tau + 0.5d; - return SelectInplace(data, (int)Math.Ceiling(h - 0.5d) - 1); - } + { + double h = data.Length*tau + 0.5d; + return SelectInplace(data, (int)Math.Ceiling(h - 0.5d) - 1); + } case QuantileDefinition.R2: - { - double h = data.Length * tau + 0.5d; - return (SelectInplace(data, (int) Math.Ceiling(h - 0.5d) - 1) + SelectInplace(data, (int) (h + 0.5d) - 1))*0.5d; - } + { + double h = data.Length*tau + 0.5d; + return (SelectInplace(data, (int)Math.Ceiling(h - 0.5d) - 1) + SelectInplace(data, (int)(h + 0.5d) - 1))*0.5d; + } case QuantileDefinition.R3: - { - double h = data.Length * tau; - return SelectInplace(data, (int)Math.Round(h) - 1); - } + { + double h = data.Length*tau; + return SelectInplace(data, (int)Math.Round(h) - 1); + } case QuantileDefinition.R4: - { - double h = data.Length * tau; - var hf = (int)h; - var lower = SelectInplace(data, hf - 1); - var upper = SelectInplace(data, hf); - return lower + (h - hf) * (upper - lower); - } + { + double h = data.Length*tau; + var hf = (int)h; + var lower = SelectInplace(data, hf - 1); + var upper = SelectInplace(data, hf); + return lower + (h - hf)*(upper - lower); + } case QuantileDefinition.R5: - { - double h = data.Length * tau + 0.5d; - var hf = (int)h; - var lower = SelectInplace(data, hf - 1); - var upper = SelectInplace(data, hf); - return lower + (h - hf) * (upper - lower); - } + { + double h = data.Length*tau + 0.5d; + var hf = (int)h; + var lower = SelectInplace(data, hf - 1); + var upper = SelectInplace(data, hf); + return lower + (h - hf)*(upper - lower); + } case QuantileDefinition.R6: - { - double h = (data.Length + 1) * tau; - var hf = (int)h; - var lower = SelectInplace(data, hf - 1); - var upper = SelectInplace(data, hf); - return lower + (h - hf) * (upper - lower); - } + { + double h = (data.Length + 1)*tau; + var hf = (int)h; + var lower = SelectInplace(data, hf - 1); + var upper = SelectInplace(data, hf); + return lower + (h - hf)*(upper - lower); + } case QuantileDefinition.R7: - { - double h = (data.Length - 1) * tau + 1d; - var hf = (int)h; - var lower = SelectInplace(data, hf - 1); - var upper = SelectInplace(data, hf); - return lower + (h - hf) * (upper - lower); - } + { + double h = (data.Length - 1)*tau + 1d; + var hf = (int)h; + var lower = SelectInplace(data, hf - 1); + var upper = SelectInplace(data, hf); + return lower + (h - hf)*(upper - lower); + } case QuantileDefinition.R8: - { - double h = (data.Length + 1 / 3d) * tau + 1 / 3d; - var hf = (int)h; - var lower = SelectInplace(data, hf - 1); - var upper = SelectInplace(data, hf); - return lower + (h - hf) * (upper - lower); - } + { + double h = (data.Length + 1/3d)*tau + 1/3d; + var hf = (int)h; + var lower = SelectInplace(data, hf - 1); + var upper = SelectInplace(data, hf); + return lower + (h - hf)*(upper - lower); + } case QuantileDefinition.R9: - { - double h = (data.Length + 0.25d) * tau + 0.375d; - var hf = (int)h; - var lower = SelectInplace(data, hf - 1); - var upper = SelectInplace(data, hf); - return lower + (h - hf) * (upper - lower); - } + { + double h = (data.Length + 0.25d)*tau + 0.375d; + var hf = (int)h; + var lower = SelectInplace(data, hf - 1); + var upper = SelectInplace(data, hf); + return lower + (h - hf)*(upper - lower); + } default: throw new NotSupportedException(); } @@ -579,5 +579,87 @@ namespace MathNet.Numerics.Statistics if (end <= rank) low = begin; } } + + /// + /// Evaluates the rank of each entry of the unsorted data array. + /// The rank definition can be specificed to be compatible + /// with an existing system. + /// WARNING: Works inplace and can thus causes the data array to be reordered. + /// + public static double[] RanksInplace(double[] data, RankDefinition definition = RankDefinition.Default) + { + var ranks = new double[data.Length]; + var index = new int[data.Length]; + for (int i = 0; i < index.Length; i++) + { + index[i] = i; + } + + if (definition == RankDefinition.First) + { + Sorting.SortAll(data, index); + for (int i = 0; i < ranks.Length; i++) + { + ranks[index[i]] = i + 1; + } + return ranks; + } + + Sorting.Sort(data, index); + int previousIndex = 0; + for (int i = 1; i < data.Length; i++) + { + if (Math.Abs(data[i] - data[previousIndex]) <= 0d) + { + continue; + } + + if (i == previousIndex + 1) + { + ranks[index[previousIndex]] = i; + } + else + { + RanksTies(ranks, index, previousIndex, i, definition); + } + + previousIndex = i; + } + + RanksTies(ranks, index, previousIndex, data.Length, definition); + return ranks; + } + + static void RanksTies(double[] ranks, int[] index, int a, int b, RankDefinition definition) + { + // TODO: potential for PERF optimization + + double rank; + switch (definition) + { + case RankDefinition.Average: + { + rank = (b + a - 1)/2d + 1; + break; + } + case RankDefinition.Min: + { + rank = a + 1; + break; + } + case RankDefinition.Max: + { + rank = b; + break; + } + default: + throw new NotSupportedException(); + } + + for (int k = a; k < b; k++) + { + ranks[index[k]] = rank; + } + } } -} \ No newline at end of file +} diff --git a/src/Numerics/Statistics/Correlation.cs b/src/Numerics/Statistics/Correlation.cs index c331b34c..95b875bd 100644 --- a/src/Numerics/Statistics/Correlation.cs +++ b/src/Numerics/Statistics/Correlation.cs @@ -161,41 +161,10 @@ namespace MathNet.Numerics.Statistics } // WARNING: do not try to cast series to an array and use it directly, - // as we need to sort it (and thus modify id) + // as we need to sort it (inplace operation) - double[] samples = series.ToArray(); - int[] index = new int[samples.Length]; - for (int i = 0; i < index.Length; i++) - { - index[i] = i; - } - Sorting.Sort(samples, index); - - double[] rankedArray = new double[samples.Length]; - int previousIndex = 0; - for (int i = 1; i < samples.Length; i++) - { - if (Math.Abs(samples[i] - samples[previousIndex]) <= 0d) - { - continue; - } - - var rankedValue = (i + previousIndex - 1) / 2d + 1; - for (int k = previousIndex; k < i; k++) - { - rankedArray[index[k]] = rankedValue; - } - - previousIndex = i; - } - - var finalValue = (samples.Length + previousIndex - 1) / 2d + 1; - for (int k = previousIndex; k < index.Length; k++) - { - rankedArray[index[k]] = finalValue; - } - - return rankedArray; + var data = series.ToArray(); + return ArrayStatistics.RanksInplace(data, RankDefinition.Average); } } } diff --git a/src/Numerics/Statistics/RankDefinition.cs b/src/Numerics/Statistics/RankDefinition.cs new file mode 100644 index 00000000..a34e3825 --- /dev/null +++ b/src/Numerics/Statistics/RankDefinition.cs @@ -0,0 +1,49 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2013 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +namespace MathNet.Numerics.Statistics +{ + public enum RankDefinition + { + /// Replace ties with their mean (non-integer ranks). Default. + Average = 1, + Default = 1, + + /// Replace ties with their minimum (typical sports ranking). + Min = 2, + Sports = 2, + + /// Replace ties with their maximum. + Max = 3, + + /// Permutation with increasing values at each index of ties. + First = 4 + } +} diff --git a/src/Numerics/Statistics/SortedArrayStatistics.cs b/src/Numerics/Statistics/SortedArrayStatistics.cs index 347b97cd..95397227 100644 --- a/src/Numerics/Statistics/SortedArrayStatistics.cs +++ b/src/Numerics/Statistics/SortedArrayStatistics.cs @@ -133,8 +133,8 @@ namespace MathNet.Numerics.Statistics /// Sample array, must be sorted ascendingly. public static double[] FiveNumberSummary(double[] data) { - if (data.Length == 0) return new[] {double.NaN, double.NaN, double.NaN, double.NaN, double.NaN}; - return new[] {data[0], Quantile(data, 0.25), Quantile(data, 0.50), Quantile(data, 0.75), data[data.Length - 1]}; + if (data.Length == 0) return new[] { double.NaN, double.NaN, double.NaN, double.NaN, double.NaN }; + return new[] { data[0], Quantile(data, 0.25), Quantile(data, 0.50), Quantile(data, 0.75), data[data.Length - 1] }; } /// @@ -157,10 +157,10 @@ namespace MathNet.Numerics.Statistics if (tau == 1d) return data[data.Length - 1]; double h = (data.Length + 1/3d)*tau + 1/3d; - var hf = (int) h; + var hf = (int)h; return hf < 1 ? data[0] : hf >= data.Length ? data[data.Length - 1] - : data[hf - 1] + (h - hf)*(data[hf] - data[hf - 1]); + : data[hf - 1] + (h - hf)*(data[hf] - data[hf - 1]); } /// @@ -189,11 +189,11 @@ namespace MathNet.Numerics.Statistics if (Math.Abs(fp) < 1e-9) { - return data[Math.Min(Math.Max((int) ip, 0), data.Length - 1)]; + return data[Math.Min(Math.Max((int)ip, 0), data.Length - 1)]; } - var lower = data[Math.Max((int) Math.Floor(x), 0)]; - var upper = data[Math.Min((int) Math.Ceiling(x), data.Length - 1)]; + var lower = data[Math.Max((int)Math.Floor(x), 0)]; + var upper = data[Math.Min((int)Math.Ceiling(x), data.Length - 1)]; return lower + (upper - lower)*(c + d*fp); } @@ -215,71 +215,145 @@ namespace MathNet.Numerics.Statistics switch (definition) { case QuantileDefinition.R1: - { - double h = data.Length*tau + 0.5d; - return data[(int) Math.Ceiling(h - 0.5d) - 1]; - } + { + double h = data.Length*tau + 0.5d; + return data[(int)Math.Ceiling(h - 0.5d) - 1]; + } case QuantileDefinition.R2: - { - double h = data.Length*tau + 0.5d; - return (data[(int) Math.Ceiling(h - 0.5d) - 1] + data[(int) (h + 0.5d) - 1])*0.5d; - } + { + double h = data.Length*tau + 0.5d; + return (data[(int)Math.Ceiling(h - 0.5d) - 1] + data[(int)(h + 0.5d) - 1])*0.5d; + } case QuantileDefinition.R3: - { - double h = data.Length*tau; - return data[Math.Max((int) Math.Round(h) - 1, 0)]; - } + { + double h = data.Length*tau; + return data[Math.Max((int)Math.Round(h) - 1, 0)]; + } case QuantileDefinition.R4: - { - double h = data.Length*tau; - var hf = (int) h; - var lower = data[Math.Max(hf - 1, 0)]; - var upper = data[Math.Min(hf, data.Length - 1)]; - return lower + (h - hf)*(upper - lower); - } + { + double h = data.Length*tau; + var hf = (int)h; + var lower = data[Math.Max(hf - 1, 0)]; + var upper = data[Math.Min(hf, data.Length - 1)]; + return lower + (h - hf)*(upper - lower); + } case QuantileDefinition.R5: - { - double h = data.Length*tau + 0.5d; - var hf = (int) h; - var lower = data[Math.Max(hf - 1, 0)]; - var upper = data[Math.Min(hf, data.Length - 1)]; - return lower + (h - hf)*(upper - lower); - } + { + double h = data.Length*tau + 0.5d; + var hf = (int)h; + var lower = data[Math.Max(hf - 1, 0)]; + var upper = data[Math.Min(hf, data.Length - 1)]; + return lower + (h - hf)*(upper - lower); + } case QuantileDefinition.R6: - { - double h = (data.Length + 1)*tau; - var hf = (int) h; - var lower = data[Math.Max(hf - 1, 0)]; - var upper = data[Math.Min(hf, data.Length - 1)]; - return lower + (h - hf)*(upper - lower); - } + { + double h = (data.Length + 1)*tau; + var hf = (int)h; + var lower = data[Math.Max(hf - 1, 0)]; + var upper = data[Math.Min(hf, data.Length - 1)]; + return lower + (h - hf)*(upper - lower); + } case QuantileDefinition.R7: - { - double h = (data.Length - 1)*tau + 1d; - var hf = (int) h; - var lower = data[Math.Max(hf - 1, 0)]; - var upper = data[Math.Min(hf, data.Length - 1)]; - return lower + (h - hf)*(upper - lower); - } + { + double h = (data.Length - 1)*tau + 1d; + var hf = (int)h; + var lower = data[Math.Max(hf - 1, 0)]; + var upper = data[Math.Min(hf, data.Length - 1)]; + return lower + (h - hf)*(upper - lower); + } case QuantileDefinition.R8: - { - double h = (data.Length + 1/3d)*tau + 1/3d; - var hf = (int) h; - var lower = data[Math.Max(hf - 1, 0)]; - var upper = data[Math.Min(hf, data.Length - 1)]; - return lower + (h - hf)*(upper - lower); - } + { + double h = (data.Length + 1/3d)*tau + 1/3d; + var hf = (int)h; + var lower = data[Math.Max(hf - 1, 0)]; + var upper = data[Math.Min(hf, data.Length - 1)]; + return lower + (h - hf)*(upper - lower); + } case QuantileDefinition.R9: - { - double h = (data.Length + 0.25d)*tau + 0.375d; - var hf = (int) h; - var lower = data[Math.Max(hf - 1, 0)]; - var upper = data[Math.Min(hf, data.Length - 1)]; - return lower + (h - hf)*(upper - lower); - } + { + double h = (data.Length + 0.25d)*tau + 0.375d; + var hf = (int)h; + var lower = data[Math.Max(hf - 1, 0)]; + var upper = data[Math.Min(hf, data.Length - 1)]; + return lower + (h - hf)*(upper - lower); + } default: throw new NotSupportedException(); } } + + /// + /// Evaluates the rank of each entry of the sorted data array (ascending). + /// The rank definition can be specificed to be compatible + /// with an existing system. + /// + public static double[] Ranks(double[] data, RankDefinition definition = RankDefinition.Default) + { + var ranks = new double[data.Length]; + + if (definition == RankDefinition.First) + { + for (int i = 0; i < ranks.Length; i++) + { + ranks[i] = i + 1; + } + return ranks; + } + + int previousIndex = 0; + for (int i = 1; i < data.Length; i++) + { + if (Math.Abs(data[i] - data[previousIndex]) <= 0d) + { + continue; + } + + if (i == previousIndex + 1) + { + ranks[previousIndex] = i; + } + else + { + RanksTies(ranks, previousIndex, i, definition); + } + + previousIndex = i; + } + + RanksTies(ranks, previousIndex, data.Length, definition); + return ranks; + } + + static void RanksTies(double[] ranks, int a, int b, RankDefinition definition) + { + // TODO: potential for PERF optimization + + double rank; + switch (definition) + { + case RankDefinition.Average: + { + rank = (b + a - 1)/2d + 1; + break; + } + case RankDefinition.Min: + { + rank = a + 1; + break; + } + case RankDefinition.Max: + { + rank = b; + break; + } + default: + throw new NotSupportedException(); + } + + for (int k = a; k < b; k++) + { + ranks[k] = rank; + } + } } } diff --git a/src/Numerics/Statistics/Statistics.cs b/src/Numerics/Statistics/Statistics.cs index 73369c73..b10c0c13 100644 --- a/src/Numerics/Statistics/Statistics.cs +++ b/src/Numerics/Statistics/Statistics.cs @@ -634,5 +634,31 @@ namespace MathNet.Numerics.Statistics Array.Sort(array); return order => SortedArrayStatistics.OrderStatistic(array, order); } + + /// + /// Evaluates the rank of each entry of the provided samples. + /// The rank definition can be specificed to be compatible + /// with an existing system. + /// + /// The data sample sequence. + /// Rank definition, to choose how ties should be handled. + public static double[] Ranks(this IEnumerable data, RankDefinition definition = RankDefinition.Default) + { + var array = data.ToArray(); + return ArrayStatistics.RanksInplace(array, definition); + } + + /// + /// Evaluates the rank of each entry of the provided samples. + /// The rank definition can be specificed to be compatible + /// with an existing system. + /// + /// The data sample sequence. + /// Rank definition, to choose how ties should be handled. + public static double[] Ranks(this IEnumerable data, RankDefinition definition = RankDefinition.Default) + { + var array = data.Where(d => d.HasValue).Select(d => d.Value).ToArray(); + return ArrayStatistics.RanksInplace(array, definition); + } } } diff --git a/src/UnitTests/ExcelTests.cs b/src/UnitTests/ExcelTests.cs index e34e41ec..1280a5b8 100644 --- a/src/UnitTests/ExcelTests.cs +++ b/src/UnitTests/ExcelTests.cs @@ -94,6 +94,13 @@ namespace MathNet.Numerics.UnitTests Assert.That(ExcelFunctions.QUARTILE(array, 2), Is.EqualTo(7.50000000000).Within(1e-8)); Assert.That(ExcelFunctions.QUARTILE(array, 3), Is.EqualTo(9.25000000000).Within(1e-8)); Assert.That(ExcelFunctions.QUARTILE(array, 4), Is.EqualTo(12.00000000000).Within(1e-8)); + + array = new Double[] { 1, 9, 12, 7, 2, 9, 10, 2 }; + Assert.That(ExcelFunctions.QUARTILE(array, 0), Is.EqualTo(1.00000000000).Within(1e-8)); + Assert.That(ExcelFunctions.QUARTILE(array, 1), Is.EqualTo(2.00000000000).Within(1e-8)); + Assert.That(ExcelFunctions.QUARTILE(array, 2), Is.EqualTo(8.00000000000).Within(1e-8)); + Assert.That(ExcelFunctions.QUARTILE(array, 3), Is.EqualTo(9.25000000000).Within(1e-8)); + Assert.That(ExcelFunctions.QUARTILE(array, 4), Is.EqualTo(12.00000000000).Within(1e-8)); } } } diff --git a/src/UnitTests/StatisticsTests/StatisticsTests.cs b/src/UnitTests/StatisticsTests/StatisticsTests.cs index b0dd0a38..725d8f4e 100644 --- a/src/UnitTests/StatisticsTests/StatisticsTests.cs +++ b/src/UnitTests/StatisticsTests/StatisticsTests.cs @@ -36,6 +36,8 @@ using MathNet.Numerics.Distributions; using MathNet.Numerics.Random; using NUnit.Framework; +// ReSharper disable InvokeAsExtensionMethod + namespace MathNet.Numerics.UnitTests.StatisticsTests { using Statistics; @@ -540,6 +542,123 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests Assert.AreEqual(expected, SortedArrayStatistics.QuantileCustom(samples, tau, 3/8d, 1/4d, 0d, 1d), 1e-14); } + [Test] + public void RanksSortedArray() + { + var distinct = new double[] { 1, 2, 4, 7, 8, 9, 10, 12 }; + var ties = new double[] { 1, 2, 2, 7, 9, 9, 10, 12 }; + + // R: rank(sort(data), ties.method="average") + Assert.That( + SortedArrayStatistics.Ranks(distinct, RankDefinition.Average), + Is.EqualTo(new[] { 1.0, 2, 3, 4, 5, 6, 7, 8 }).AsCollection.Within(1e-8)); + Assert.That( + SortedArrayStatistics.Ranks(ties, RankDefinition.Average), + Is.EqualTo(new[] { 1, 2.5, 2.5, 4, 5.5, 5.5, 7, 8 }).AsCollection.Within(1e-8)); + + // R: rank(data, ties.method="min") + Assert.That( + SortedArrayStatistics.Ranks(distinct, RankDefinition.Min), + Is.EqualTo(new[] { 1.0, 2, 3, 4, 5, 6, 7, 8 }).AsCollection.Within(1e-8)); + Assert.That( + SortedArrayStatistics.Ranks(ties, RankDefinition.Min), + Is.EqualTo(new[] { 1.0, 2, 2, 4, 5, 5, 7, 8 }).AsCollection.Within(1e-8)); + + // R: rank(data, ties.method="max") + Assert.That( + SortedArrayStatistics.Ranks(distinct, RankDefinition.Max), + Is.EqualTo(new[] { 1.0, 2, 3, 4, 5, 6, 7, 8 }).AsCollection.Within(1e-8)); + Assert.That( + SortedArrayStatistics.Ranks(ties, RankDefinition.Max), + Is.EqualTo(new[] { 1.0, 3, 3, 4, 6, 6, 7, 8 }).AsCollection.Within(1e-8)); + + // R: rank(data, ties.method="first") + Assert.That( + SortedArrayStatistics.Ranks(distinct, RankDefinition.First), + Is.EqualTo(new[] { 1.0, 2, 3, 4, 5, 6, 7, 8 }).AsCollection.Within(1e-8)); + Assert.That( + SortedArrayStatistics.Ranks(ties, RankDefinition.First), + Is.EqualTo(new[] { 1.0, 2, 3, 4, 5, 6, 7, 8 }).AsCollection.Within(1e-8)); + } + + [Test] + public void RanksArray() + { + var distinct = new double[] { 1, 8, 12, 7, 2, 9, 10, 4 }; + var ties = new double[] { 1, 9, 12, 7, 2, 9, 10, 2 }; + + // R: rank(data, ties.method="average") + Assert.That( + ArrayStatistics.RanksInplace((double[])distinct.Clone(), RankDefinition.Average), + Is.EqualTo(new[] { 1.0, 5, 8, 4, 2, 6, 7, 3 }).AsCollection.Within(1e-8)); + Assert.That( + ArrayStatistics.RanksInplace((double[])ties.Clone(), RankDefinition.Average), + Is.EqualTo(new[] { 1, 5.5, 8, 4, 2.5, 5.5, 7, 2.5 }).AsCollection.Within(1e-8)); + + // R: rank(data, ties.method="min") + Assert.That( + ArrayStatistics.RanksInplace((double[])distinct.Clone(), RankDefinition.Min), + Is.EqualTo(new[] { 1.0, 5, 8, 4, 2, 6, 7, 3 }).AsCollection.Within(1e-8)); + Assert.That( + ArrayStatistics.RanksInplace((double[])ties.Clone(), RankDefinition.Min), + Is.EqualTo(new[] { 1.0, 5, 8, 4, 2, 5, 7, 2 }).AsCollection.Within(1e-8)); + + // R: rank(data, ties.method="max") + Assert.That( + ArrayStatistics.RanksInplace((double[])distinct.Clone(), RankDefinition.Max), + Is.EqualTo(new[] { 1.0, 5, 8, 4, 2, 6, 7, 3 }).AsCollection.Within(1e-8)); + Assert.That( + ArrayStatistics.RanksInplace((double[])ties.Clone(), RankDefinition.Max), + Is.EqualTo(new[] { 1.0, 6, 8, 4, 3, 6, 7, 3 }).AsCollection.Within(1e-8)); + + // R: rank(data, ties.method="first") + Assert.That( + ArrayStatistics.RanksInplace((double[])distinct.Clone(), RankDefinition.First), + Is.EqualTo(new[] { 1.0, 5, 8, 4, 2, 6, 7, 3 }).AsCollection.Within(1e-8)); + Assert.That( + ArrayStatistics.RanksInplace((double[])ties.Clone(), RankDefinition.First), + Is.EqualTo(new[] { 1.0, 5, 8, 4, 2, 6, 7, 3 }).AsCollection.Within(1e-8)); + } + + [Test] + public void Ranks() + { + var distinct = new double[] { 1, 8, 12, 7, 2, 9, 10, 4 }; + var ties = new double[] { 1, 9, 12, 7, 2, 9, 10, 2 }; + + // R: rank(data, ties.method="average") + Assert.That( + Statistics.Ranks(distinct, RankDefinition.Average), + Is.EqualTo(new[] { 1.0, 5, 8, 4, 2, 6, 7, 3 }).AsCollection.Within(1e-8)); + Assert.That( + Statistics.Ranks(ties, RankDefinition.Average), + Is.EqualTo(new[] { 1, 5.5, 8, 4, 2.5, 5.5, 7, 2.5 }).AsCollection.Within(1e-8)); + + // R: rank(data, ties.method="min") + Assert.That( + Statistics.Ranks(distinct, RankDefinition.Min), + Is.EqualTo(new[] { 1.0, 5, 8, 4, 2, 6, 7, 3 }).AsCollection.Within(1e-8)); + Assert.That( + Statistics.Ranks(ties, RankDefinition.Min), + Is.EqualTo(new[] { 1.0, 5, 8, 4, 2, 5, 7, 2 }).AsCollection.Within(1e-8)); + + // R: rank(data, ties.method="max") + Assert.That( + Statistics.Ranks(distinct, RankDefinition.Max), + Is.EqualTo(new[] { 1.0, 5, 8, 4, 2, 6, 7, 3 }).AsCollection.Within(1e-8)); + Assert.That( + Statistics.Ranks(ties, RankDefinition.Max), + Is.EqualTo(new[] { 1.0, 6, 8, 4, 3, 6, 7, 3 }).AsCollection.Within(1e-8)); + + // R: rank(data, ties.method="first") + Assert.That( + Statistics.Ranks(distinct, RankDefinition.First), + Is.EqualTo(new[] { 1.0, 5, 8, 4, 2, 6, 7, 3 }).AsCollection.Within(1e-8)); + Assert.That( + Statistics.Ranks(ties, RankDefinition.First), + Is.EqualTo(new[] { 1.0, 5, 8, 4, 2, 6, 7, 3 }).AsCollection.Within(1e-8)); + } + [Test] public void MedianOnShortSequence() { @@ -734,3 +853,5 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests } } } + +// ReSharper restore InvokeAsExtensionMethod