From 6647d007ab896ca6e141de4361b2ab36509e49e3 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Fri, 22 Mar 2013 23:16:03 +0100 Subject: [PATCH] Statistics: Proper Quantile statistics (supporting all 9 R-types and custom) --- src/Numerics/Numerics.csproj | 1 + src/Numerics/Statistics/ArrayStatistics.cs | 263 +++++++++++++++ src/Numerics/Statistics/Histogram.cs | 9 +- src/Numerics/Statistics/Percentile.cs | 8 +- src/Numerics/Statistics/QuantileDefinition.cs | 45 +++ .../Statistics/SortedArrayStatistics.cs | 144 ++++++--- src/Numerics/Statistics/Statistics.cs | 231 ++++++------- .../Statistics/StreamingStatistics.cs | 7 + src/Portable/Portable.csproj | 3 + .../StatisticsTests/StatisticsTests.cs | 306 +++++++++++++++++- 10 files changed, 819 insertions(+), 198 deletions(-) create mode 100644 src/Numerics/Statistics/QuantileDefinition.cs diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index c3f0fd99..60be540e 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -112,6 +112,7 @@ + diff --git a/src/Numerics/Statistics/ArrayStatistics.cs b/src/Numerics/Statistics/ArrayStatistics.cs index 6e41011d..628e53ec 100644 --- a/src/Numerics/Statistics/ArrayStatistics.cs +++ b/src/Numerics/Statistics/ArrayStatistics.cs @@ -32,6 +32,13 @@ using System; namespace MathNet.Numerics.Statistics { + /// + /// Statistics operating on arrays assumed to be unsorted. + /// WARNING: Methods with the Inplace-suffix may modify the data array by reordering its entries. + /// + /// + /// + /// public static class ArrayStatistics { // TODO: Benchmark various options to find out the best approach (-> branch prediction) @@ -79,6 +86,22 @@ namespace MathNet.Numerics.Statistics return max; } + /// + /// Returns the order statistic (order 1..N) from the unsorted data array. + /// WARNING: Works inplace and can thus causes the data array to be reordered. + /// + /// Sample array, no sorting is assumed. Will be reordered. + /// One-based order of the statistic, must be between 1 and N (inclusive). + public static double OrderStatisticInplace(double[] data, int order) + { + if (data == null) throw new ArgumentNullException("data"); + if (order < 1 || order > data.Length) return double.NaN; + + if (order == 1) return Minimum(data); + if (order == data.Length) return Maximum(data); + return SelectInplace(data, order - 1); + } + /// /// Estimates the arithmetic sample mean from the unsorted data array. /// Returns NaN if data is empty or any entry is NaN. @@ -163,5 +186,245 @@ namespace MathNet.Numerics.Statistics { return Math.Sqrt(PopulationVariance(data)); } + + /// + /// Estimates the median value from the unsorted data array. + /// Applies a linear interpolation, consistent with Quantile and R-8. + /// WARNING: Works inplace and can thus causes the data array to be reordered. + /// + /// Sample array, no sorting is assumed. Will be reordered. + public static double MedianInplace(double[] data) + { + return QuantileInplace(data, 0.5d); + } + + /// + /// Estimates the tau-th quantile from the unsorted data array. + /// The tau-th quantile is the data value where the cumulative distribution + /// function crosses tau. Applies a linear interpolation, compatible with R-8. + /// WARNING: Works inplace and can thus causes the data array to be reordered. + /// + /// Sample array, no sorting is assumed. Will be reordered. + /// Quantile selector, between 0.0 and 1.0 (inclusive). + /// + /// R-8, SciPy-(1/3,1/3): + /// Linear interpolation of the approximate medians for order statistics. + /// When tau < (2/3) / (N + 1/3), use x1. When tau >= (N - 1/3) / (N + 1/3), use xN. + /// + public static double QuantileInplace(double[] data, double tau) + { + if (data == null) throw new ArgumentNullException("data"); + if (tau < 0d || tau > 1d || data.Length == 0) return double.NaN; + + double h = (data.Length + 1d/3d)*tau + 1d/3d; + var hf = (int) h; + + if (hf <= 0 || tau == 0d) + { + return Minimum(data); + } + if (hf >= data.Length || tau == 1d) + { + return Maximum(data); + } + + var a = SelectInplace(data, hf - 1); + var b = SelectInplace(data, hf); + return a + (h - hf)*(b - a); + } + + /// + /// Estimates the tau-th quantile from the unsorted data array. + /// The tau-th quantile is the data value where the cumulative distribution + /// function crosses tau. The quantile defintion can be specified + /// by 4 parameters a, b, c and d, consistent with Mathematica. + /// WARNING: Works inplace and can thus causes the data array to be reordered. + /// + /// Sample array, no sorting is assumed. Will be reordered. + /// Quantile selector, between 0.0 and 1.0 (inclusive) + public static double QuantileCustomInplace(double[] data, double tau, double a, double b, double c, double d) + { + if (data == null) throw new ArgumentNullException("data"); + if (tau < 0d || tau > 1d || data.Length == 0) return double.NaN; + + var x = a + (data.Length + b) * tau - 1; +#if PORTABLE + var ip = (int)x; +#else + var ip = Math.Truncate(x); +#endif + var fp = x - ip; + + if (Math.Abs(fp) < 1e-9) + { + 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); + } + + /// + /// Estimates the tau-th quantile from the unsorted data array. + /// The tau-th quantile is the data value where the cumulative distribution + /// function crosses tau. The quantile definition can be specificed to be compatible + /// with an existing system. + /// WARNING: Works inplace and can thus causes the data array to be reordered. + /// + /// Sample array, no sorting is assumed. Will be reordered. + /// Quantile selector, between 0.0 and 1.0 (inclusive) + /// Quantile definition, to choose what product/definition it should be consistent with + public static double QuantileCustomInplace(double[] data, double tau, QuantileDefinition definition) + { + if (data == null) throw new ArgumentNullException("data"); + if (tau < 0d || tau > 1d || data.Length == 0) return double.NaN; + if (tau == 0d || data.Length == 1) return Minimum(data); + if (tau == 1d) return Maximum(data); + + switch (definition) + { + case QuantileDefinition.R1: + { + 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; + } + case QuantileDefinition.R3: + { + 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); + } + 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); + } + 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); + } + 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); + } + 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); + } + 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); + } + default: + throw new NotSupportedException(); + } + } + + static double SelectInplace(double[] workingData, int rank) + { + // Numerical Recipes: select + // http://en.wikipedia.org/wiki/Selection_algorithm + + if (rank <= 0) return Minimum(workingData); + if (rank >= workingData.Length - 1) return Maximum(workingData); + + var a = workingData; + int low = 0; + int high = a.Length - 1; + + while (true) + { + if (high <= low + 1) + { + if (high == low + 1 && a[high] < a[low]) + { + var tmp = a[low]; + a[low] = a[high]; + a[high] = tmp; + } + + return a[rank]; + } + + int middle = (low + high) >> 1; + + var tmp1 = a[middle]; + a[middle] = a[low + 1]; + a[low + 1] = tmp1; + + if (a[low] > a[high]) + { + var tmp = a[low]; + a[low] = a[high]; + a[high] = tmp; + } + if (a[low + 1] > a[high]) + { + var tmp = a[low + 1]; + a[low + 1] = a[high]; + a[high] = tmp; + } + if (a[low] > a[low + 1]) + { + var tmp = a[low]; + a[low] = a[low + 1]; + a[low + 1] = tmp; + } + + int begin = low + 1; + int end = high; + double pivot = a[begin]; + + while (true) + { + do begin++; while (a[begin] < pivot); + do end--; while (a[end] > pivot); + + if (end < begin) break; + + var tmp = a[begin]; + a[begin] = a[end]; + a[end] = tmp; + } + + a[low + 1] = a[end]; + a[end] = pivot; + + if (end >= rank) high = end - 1; + if (end <= rank) low = begin; + } + } } } \ No newline at end of file diff --git a/src/Numerics/Statistics/Histogram.cs b/src/Numerics/Statistics/Histogram.cs index 124a774a..86aa88f2 100644 --- a/src/Numerics/Statistics/Histogram.cs +++ b/src/Numerics/Statistics/Histogram.cs @@ -246,7 +246,7 @@ namespace MathNet.Numerics.Statistics /// /// The datasequence to build a histogram on. /// The number of buckets to use. - public Histogram(IEnumerable data, int nbuckets) + public Histogram(IEnumerable data, int nbuckets) : this() { if (nbuckets < 1) @@ -256,7 +256,12 @@ namespace MathNet.Numerics.Statistics double lower = data.Minimum(); double upper = data.Maximum(); - double width = (upper - lower) / nbuckets; + double width = (upper - lower)/nbuckets; + + if (double.IsNaN(width)) + { + throw new ArgumentException("Data must contain at least one entry.", "data"); + } // Add buckets for each bin; the smallest bucket's lowerbound must be slightly smaller // than the minimal element. diff --git a/src/Numerics/Statistics/Percentile.cs b/src/Numerics/Statistics/Percentile.cs index 65f3bb8a..59c8aa29 100644 --- a/src/Numerics/Statistics/Percentile.cs +++ b/src/Numerics/Statistics/Percentile.cs @@ -109,13 +109,13 @@ namespace MathNet.Numerics.Statistics switch (Method) { case PercentileMethod.Nist: - return SortedArrayStatistics.QuantileCompatible(_data, percentile, QuantileCompatibility.Nist); + return SortedArrayStatistics.QuantileCustom(_data, percentile, QuantileDefinition.Nist); case PercentileMethod.Nearest: - return SortedArrayStatistics.QuantileCompatible(_data, percentile, QuantileCompatibility.R3); + return SortedArrayStatistics.QuantileCustom(_data, percentile, QuantileDefinition.R3); case PercentileMethod.Interpolation: - return SortedArrayStatistics.QuantileCompatible(_data, percentile, QuantileCompatibility.R5); + return SortedArrayStatistics.QuantileCustom(_data, percentile, QuantileDefinition.R5); case PercentileMethod.Excel: - return SortedArrayStatistics.QuantileCompatible(_data, percentile, QuantileCompatibility.Excel); + return SortedArrayStatistics.QuantileCustom(_data, percentile, QuantileDefinition.Excel); default: return SortedArrayStatistics.Quantile(_data, percentile); } diff --git a/src/Numerics/Statistics/QuantileDefinition.cs b/src/Numerics/Statistics/QuantileDefinition.cs new file mode 100644 index 00000000..6e7deeb7 --- /dev/null +++ b/src/Numerics/Statistics/QuantileDefinition.cs @@ -0,0 +1,45 @@ +// +// 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 QuantileDefinition + { + R1 = 1, SAS3 = 1, InverseCDF = 1, + R2 = 2, SAS5 = 2, InverseCDFAverage = 2, + R3 = 3, SAS2 = 3, Nearest = 3, + R4 = 4, SAS1 = 4, California = 4, + R5 = 5, Hydrology = 5, Hazen = 5, + R6 = 6, SAS4 = 6, Nist = 6, Weibull = 6, SPSS = 6, + R7 = 7, Excel = 7, Mode = 7, S = 7, + R8 = 8, Median = 8, Default = 8, + R9 = 9, Normal = 9, + } +} diff --git a/src/Numerics/Statistics/SortedArrayStatistics.cs b/src/Numerics/Statistics/SortedArrayStatistics.cs index fe9ab77c..7061ab8e 100644 --- a/src/Numerics/Statistics/SortedArrayStatistics.cs +++ b/src/Numerics/Statistics/SortedArrayStatistics.cs @@ -32,19 +32,14 @@ using System; namespace MathNet.Numerics.Statistics { - public enum QuantileCompatibility - { - Default=0, - Nist,Nearest,Excel, - R1,R2,R3,R4,R5,R6,R7,R8,R9, - SAS1,SAS2,SAS3,SAS4,SAS5 - } - + /// + /// Statistics operating on an array already sorted ascendingly. + /// + /// + /// + /// public static class SortedArrayStatistics { - const double Third = 1d/3d; - const double Half = 1d/2d; - /// /// Returns the smallest value from the sorted data array (ascending). /// @@ -69,6 +64,19 @@ namespace MathNet.Numerics.Statistics return data[data.Length - 1]; } + /// + /// Returns the order statistic (order 1..N) from the sorted data array (ascending). + /// + /// Sample array, must be sorted ascendingly. + /// One-based order of the statistic, must be between 1 and N (inclusive). + public static double OrderStatistic(double[] data, int order) + { + if (data == null) throw new ArgumentNullException("data"); + if (order < 1 || order > data.Length) return double.NaN; + + return data[order - 1]; + } + /// /// Estimates the median value from the sorted data array (ascending). /// Applies a linear interpolation, consistent with Quantile and R-8. @@ -152,84 +160,124 @@ namespace MathNet.Numerics.Statistics if (tau == 0d || data.Length == 1) return data[0]; if (tau == 1d) return data[data.Length - 1]; - double h = (data.Length + Third)*tau + Third; + double h = (data.Length + 1/3d)*tau + 1/3d; var hf = (int) h; - return data[hf - 1] + (h - hf)*(data[hf] - data[hf - 1]); + return hf < 1 ? data[0] + : hf >= data.Length ? data[data.Length - 1] + : data[hf - 1] + (h - hf)*(data[hf] - data[hf - 1]); + } + + /// + /// Estimates the tau-th quantile from the sorted data array (ascending). + /// The tau-th quantile is the data value where the cumulative distribution + /// function crosses tau. The quantile defintion can be specified + /// by 4 parameters a, b, c and d, consistent with Mathematica. + /// + /// Sample array, must be sorted ascendingly. + /// Quantile selector, between 0.0 and 1.0 (inclusive). + public static double QuantileCustom(double[] data, double tau, double a, double b, double c, double d) + { + if (data == null) throw new ArgumentNullException("data"); + if (tau < 0d || tau > 1d || data.Length == 0) return double.NaN; + + var x = a + (data.Length + b)*tau - 1; +#if PORTABLE + var ip = (int) x; +#else + var ip = Math.Truncate(x); +#endif + var fp = x - ip; + + if (Math.Abs(fp) < 1e-9) + { + 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)]; + return lower + (upper - lower)*(c + d*fp); } /// /// Estimates the tau-th quantile from the sorted data array (ascending). /// The tau-th quantile is the data value where the cumulative distribution - /// function crosses tau. The quantile algorithm can be chosen by the compatibility argument. + /// function crosses tau. The quantile definition can be specificed to be compatible + /// with an existing system. /// - public static double QuantileCompatible(double[] data, double tau, QuantileCompatibility compatibility) + /// Sample array, must be sorted ascendingly. + /// Quantile selector, between 0.0 and 1.0 (inclusive). + /// Quantile definition, to choose what product/definition it should be consistent with + public static double QuantileCustom(double[] data, double tau, QuantileDefinition definition) { if (data == null) throw new ArgumentNullException("data"); if (tau < 0d || tau > 1d || data.Length == 0) return double.NaN; if (tau == 0d || data.Length == 1) return data[0]; if (tau == 1d) return data[data.Length - 1]; - switch (compatibility) + switch (definition) { - case QuantileCompatibility.R1: - case QuantileCompatibility.SAS3: + case QuantileDefinition.R1: { - double h = data.Length*tau + Half; - return data[(int) Math.Ceiling(h - Half) - 1]; + double h = data.Length*tau + 0.5d; + return data[(int) Math.Ceiling(h - 0.5d) - 1]; } - case QuantileCompatibility.R2: - case QuantileCompatibility.SAS5: + case QuantileDefinition.R2: { - double h = data.Length*tau + Half; - return (data[(int) Math.Ceiling(h - Half) - 1] + data[(int) (h + Half) - 1])*Half; + 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 QuantileCompatibility.R3: - case QuantileCompatibility.SAS2: - case QuantileCompatibility.Nearest: + case QuantileDefinition.R3: { double h = data.Length*tau; - return data[(int) Math.Round(h) - 1]; + return data[Math.Max((int) Math.Round(h) - 1, 0)]; } - case QuantileCompatibility.R4: - case QuantileCompatibility.SAS1: + case QuantileDefinition.R4: { double h = data.Length*tau; var hf = (int) h; - return data[hf - 1] + (h - hf)*(data[hf] - data[hf - 1]); + var lower = data[Math.Max(hf - 1, 0)]; + var upper = data[Math.Min(hf, data.Length - 1)]; + return lower + (h - hf)*(upper - lower); } - case QuantileCompatibility.R5: + case QuantileDefinition.R5: { - double h = data.Length*tau + Half; + double h = data.Length*tau + 0.5d; var hf = (int) h; - return data[hf - 1] + (h - hf)*(data[hf] - data[hf - 1]); + var lower = data[Math.Max(hf - 1, 0)]; + var upper = data[Math.Min(hf, data.Length - 1)]; + return lower + (h - hf)*(upper - lower); } - case QuantileCompatibility.R6: - case QuantileCompatibility.SAS4: - case QuantileCompatibility.Nist: + case QuantileDefinition.R6: { double h = (data.Length + 1)*tau; var hf = (int) h; - return data[hf - 1] + (h - hf)*(data[hf] - data[hf - 1]); + var lower = data[Math.Max(hf - 1, 0)]; + var upper = data[Math.Min(hf, data.Length - 1)]; + return lower + (h - hf)*(upper - lower); } - case QuantileCompatibility.R7: - case QuantileCompatibility.Excel: + case QuantileDefinition.R7: { double h = (data.Length - 1)*tau + 1d; var hf = (int) h; - return data[hf - 1] + (h - hf)*(data[hf] - data[hf - 1]); + var lower = data[Math.Max(hf - 1, 0)]; + var upper = data[Math.Min(hf, data.Length - 1)]; + return lower + (h - hf)*(upper - lower); } - case QuantileCompatibility.R8: - case QuantileCompatibility.Default: + case QuantileDefinition.R8: { - double h = (data.Length + Third)*tau + Third; + double h = (data.Length + 1/3d)*tau + 1/3d; var hf = (int) h; - return data[hf - 1] + (h - hf)*(data[hf] - data[hf - 1]); + var lower = data[Math.Max(hf - 1, 0)]; + var upper = data[Math.Min(hf, data.Length - 1)]; + return lower + (h - hf)*(upper - lower); } - case QuantileCompatibility.R9: + case QuantileDefinition.R9: { - double h = (data.Length + 1d/4d)*tau + 3d/8d; + double h = (data.Length + 0.25d)*tau + 0.375d; var hf = (int) h; - return data[hf - 1] + (h - hf)*(data[hf] - data[hf - 1]); + 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(); diff --git a/src/Numerics/Statistics/Statistics.cs b/src/Numerics/Statistics/Statistics.cs index 2917ab67..f1a41fed 100644 --- a/src/Numerics/Statistics/Statistics.cs +++ b/src/Numerics/Statistics/Statistics.cs @@ -33,7 +33,6 @@ namespace MathNet.Numerics.Statistics using System; using System.Collections.Generic; using System.Linq; - using Properties; /// /// Extension methods to return basic statistics on set of data. @@ -88,7 +87,7 @@ namespace MathNet.Numerics.Statistics } /// - /// Calculates the sample mean. + /// Estimates the sample mean. /// /// The data to calculate the mean of. /// The mean of the sample. @@ -101,7 +100,7 @@ namespace MathNet.Numerics.Statistics } /// - /// Calculates the sample mean. + /// Estimates the sample mean. /// /// The data to calculate the mean of. /// The mean of the sample. @@ -112,7 +111,7 @@ namespace MathNet.Numerics.Statistics } /// - /// Calculates the unbiased population (sample) variance estimator (on a dataset of size N will use an N-1 normalizer). + /// Estimates the unbiased population (sample) variance estimator (on a dataset of size N will use an N-1 normalizer). /// /// The data to calculate the variance of. /// The unbiased population variance of the sample. @@ -125,7 +124,7 @@ namespace MathNet.Numerics.Statistics } /// - /// Computes the unbiased population (sample) variance estimator (on a dataset of size N will use an N-1 normalizer) for nullable data. + /// Estimates the unbiased population (sample) variance estimator (on a dataset of size N will use an N-1 normalizer) for nullable data. /// /// The data to calculate the variance of. /// The population variance of the sample. @@ -136,7 +135,7 @@ namespace MathNet.Numerics.Statistics } /// - /// Calculates the biased population variance estimator (on a dataset of size N will use an N normalizer). + /// Estimates the biased population variance estimator (on a dataset of size N will use an N normalizer). /// /// The data to calculate the variance of. /// The biased population variance of the sample. @@ -149,7 +148,7 @@ namespace MathNet.Numerics.Statistics } /// - /// Computes the biased population variance estimator (on a dataset of size N will use an N normalizer) for nullable data. + /// Estimates the biased population variance estimator (on a dataset of size N will use an N normalizer) for nullable data. /// /// The data to calculate the variance of. /// The population variance of the sample. @@ -160,7 +159,7 @@ namespace MathNet.Numerics.Statistics } /// - /// Calculates the unbiased sample standard deviation (on a dataset of size N will use an N-1 normalizer). + /// Estimates the unbiased sample standard deviation (on a dataset of size N will use an N-1 normalizer). /// /// The data to calculate the standard deviation of. /// The standard deviation of the sample. @@ -173,7 +172,7 @@ namespace MathNet.Numerics.Statistics } /// - /// Calculates the unbiased sample standard deviation (on a dataset of size N will use an N-1 normalizer). + /// Estimates the unbiased sample standard deviation (on a dataset of size N will use an N-1 normalizer). /// /// The data to calculate the standard deviation of. /// The standard deviation of the sample. @@ -184,7 +183,7 @@ namespace MathNet.Numerics.Statistics } /// - /// Calculates the biased sample standard deviation (on a dataset of size N will use an N normalizer). + /// Estimates the biased sample standard deviation (on a dataset of size N will use an N normalizer). /// /// The data to calculate the standard deviation of. /// The standard deviation of the sample. @@ -197,7 +196,7 @@ namespace MathNet.Numerics.Statistics } /// - /// Calculates the biased sample standard deviation (on a dataset of size N will use an N normalizer). + /// Estimates the biased sample standard deviation (on a dataset of size N will use an N normalizer). /// /// The data to calculate the standard deviation of. /// The standard deviation of the sample. @@ -208,158 +207,120 @@ namespace MathNet.Numerics.Statistics } /// - /// Calculates the sample median. + /// Estimates the sample median. /// /// The data to calculate the median of. /// The median of the sample. public static double Median(this IEnumerable data) { - if (data == null) - { - throw new ArgumentNullException("data"); - } - - var dataArray = new List(data); - if (dataArray.Count == 0) - { - return double.NaN; - } - - int index = (dataArray.Count / 2) + 1; - if (dataArray.Count % 2 == 0) - { - double lower = OrderSelect(dataArray, 0, dataArray.Count - 1, index - 1); - double upper = dataArray.Skip(index - 1).Minimum(); - return (lower + upper) / 2.0; - } - - return OrderSelect(dataArray, 0, dataArray.Count - 1, index); + if (data == null) throw new ArgumentNullException("data"); + var array = data.ToArray(); + return ArrayStatistics.MedianInplace(array); } /// - /// Calculates the sample median. + /// Estimates the sample median. /// /// The data to calculate the median of. /// The median of the sample. public static double Median(this IEnumerable data) { - if (data == null) - { - throw new ArgumentNullException("data"); - } - - var nonNull = new List(); - foreach (double? value in data) - { - if (value.HasValue) - { - nonNull.Add(value.Value); - } - } - - return nonNull.Median(); + if (data == null) throw new ArgumentNullException("data"); + var array = data.Where(d => d.HasValue).Select(d => d.Value).ToArray(); + return ArrayStatistics.MedianInplace(array); } /// - /// Evaluate the i-order (1..N) statistic of the provided samples. + /// Estimates the sample tau-quantile. /// - /// The sample data. - /// Order of the statistic to evaluate. - /// The i'th order statistic in the sample data. - public static double OrderStatistic(IEnumerable samples, int order) + /// The data to calculate the median of. + /// Quantile selector, between 0.0 and 1.0 (inclusive). + /// The median of the sample. + public static double Quantile(this IEnumerable data, double tau) { - if (order == 1) - { - // Can be done in linear time by Min() - return Minimum(samples); - } - - var list = new List(samples); - if (list.Count == 0) - { - return double.NaN; - } - - if (order < 1 || order > list.Count) - { - throw new ArgumentOutOfRangeException("order", Resources.ArgumentInIntervalXYInclusive); - } - - if (order == list.Count) - { - // Can be done in linear time by Max() - return Maximum(list); - } - - return OrderSelect(list, 0, list.Count - 1, order); + if (data == null) throw new ArgumentNullException("data"); + var array = data.ToArray(); + return ArrayStatistics.QuantileInplace(array, tau); } /// - /// Implementation of the order statistics finding algorithm based on the algorithm in - /// "Introduction to Algorithms", Cormen et al. section 7.1. + /// Estimates the sample tau-quantile. /// - /// The sample data. - /// The left bound in which to order select. - /// The right bound in which to order select. - /// The order we are trying to find. - /// The order statistic. - private static double OrderSelect(IList samples, int left, int right, int order) + /// The data to calculate the median of. + /// Quantile selector, between 0.0 and 1.0 (inclusive). + /// The median of the sample. + public static double Quantile(this IEnumerable data, double tau) { - while (true) - { - System.Diagnostics.Debug.Assert(order > 0, "Order must always be positive."); - System.Diagnostics.Debug.Assert(left >= 0 && left <= right, "Left side must always be positive and smaller than right side."); - System.Diagnostics.Debug.Assert(right < samples.Count, "Right side must always be smaller than number of elements in list."); - System.Diagnostics.Debug.Assert(right - left + 1 >= order, "Make sure there are at least order items in the segment [left, right]."); - - if (left == right) - { - return samples[left]; - } - - - // The pivot point. Choose median of left, right and center - //to be the pivot and arrange so that - //samples[left]<=samples[right]<=samples[center] - int center = (left + right) / 2; - if (samples[center] < samples[left]) - Sorting.Swap(samples, left, center); - if (samples[center] < samples[right]) - Sorting.Swap(samples, right, center); - if (samples[right] < samples[left]) - Sorting.Swap(samples, right, left); + if (data == null) throw new ArgumentNullException("data"); + var array = data.Where(d => d.HasValue).Select(d => d.Value).ToArray(); + return ArrayStatistics.QuantileInplace(array, tau); + } - double pivot = samples[right]; + /// + /// Estimates the empiric inverse CDF at tau (tau-quantile). + /// + /// The data to calculate the median of. + /// Quantile selector, between 0.0 and 1.0 (inclusive). + /// The median of the sample. + public static double InverseCDF(this IEnumerable data, double tau) + { + if (data == null) throw new ArgumentNullException("data"); + var array = data.ToArray(); + return ArrayStatistics.QuantileCustomInplace(array, tau, QuantileDefinition.InverseCDF); + } - // The partioning code. - int i = left; - for (int j = left+1; j <= right - 1; j++) - { - if (samples[j] <= pivot) - { - i++; - Sorting.Swap(samples, i, j); - } - } + /// + /// Estimates the empiric inverse CDF at tau (tau-quantile). + /// + /// The data to calculate the median of. + /// Quantile selector, between 0.0 and 1.0 (inclusive). + /// The median of the sample. + public static double InverseCDF(this IEnumerable data, double tau) + { + if (data == null) throw new ArgumentNullException("data"); + var array = data.Where(d => d.HasValue).Select(d => d.Value).ToArray(); + return ArrayStatistics.QuantileCustomInplace(array, tau, QuantileDefinition.InverseCDF); + } - Sorting.Swap(samples, i + 1, right); + /// + /// Estimates the sample tau-quantile. + /// + /// The data to calculate the median of. + /// Quantile selector, between 0.0 and 1.0 (inclusive). + /// The median of the sample. + /// Quantile definition, to choose what product/definition it should be consistent with + public static double QuantileCustom(this IEnumerable data, double tau, QuantileDefinition definition) + { + if (data == null) throw new ArgumentNullException("data"); + var array = data.ToArray(); + return ArrayStatistics.QuantileCustomInplace(array, tau, definition); + } - // Recursive order finding algorithm. - if (order == (i - left) + 2) - { - return pivot; - } + /// + /// Estimates the sample tau-quantile. + /// + /// The data to calculate the median of. + /// Quantile selector, between 0.0 and 1.0 (inclusive). + /// The median of the sample. + /// Quantile definition, to choose what product/definition it should be consistent with + public static double QuantileCustom(this IEnumerable data, double tau, QuantileDefinition definition) + { + if (data == null) throw new ArgumentNullException("data"); + var array = data.Where(d => d.HasValue).Select(d => d.Value).ToArray(); + return ArrayStatistics.QuantileCustomInplace(array, tau, definition); + } - if (order < (i - left) + 2) - { - right = i; - } - else - { - order = order - i + left - 2; - left = i + 2; - } - } + /// + /// Returns the i-order (1..N) statistic of the provided samples. + /// + /// The sample data. + /// Order of the statistic to evaluate. + /// The i'th order statistic in the sample data. + public static double OrderStatistic(IEnumerable data, int order) + { + if (data == null) throw new ArgumentNullException("data"); + var array = data.ToArray(); + return ArrayStatistics.OrderStatisticInplace(array, order); } } } diff --git a/src/Numerics/Statistics/StreamingStatistics.cs b/src/Numerics/Statistics/StreamingStatistics.cs index 6a60022f..f264be90 100644 --- a/src/Numerics/Statistics/StreamingStatistics.cs +++ b/src/Numerics/Statistics/StreamingStatistics.cs @@ -33,6 +33,13 @@ using System.Collections.Generic; namespace MathNet.Numerics.Statistics { + /// + /// Statistics operating on an IEnumerable in a single pass, without keeping the full data in memory. + /// Can be used in a streaming way, e.g. on large datasets not fitting into memory. + /// + /// + /// + /// public static class StreamingStatistics { /// diff --git a/src/Portable/Portable.csproj b/src/Portable/Portable.csproj index dbf23f19..d279b8b5 100644 --- a/src/Portable/Portable.csproj +++ b/src/Portable/Portable.csproj @@ -1062,6 +1062,9 @@ Statistics\Percentile.cs + + Statistics\QuantileDefinition.cs + Statistics\SortedArrayStatistics.cs diff --git a/src/UnitTests/StatisticsTests/StatisticsTests.cs b/src/UnitTests/StatisticsTests/StatisticsTests.cs index fe7cfcf7..f06d468d 100644 --- a/src/UnitTests/StatisticsTests/StatisticsTests.cs +++ b/src/UnitTests/StatisticsTests/StatisticsTests.cs @@ -64,6 +64,7 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests Assert.Throws(() => Statistics.Maximum(data)); Assert.Throws(() => Statistics.Mean(data)); Assert.Throws(() => Statistics.Median(data)); + Assert.Throws(() => Statistics.Quantile(data, 0.3)); Assert.Throws(() => Statistics.Variance(data)); Assert.Throws(() => Statistics.StandardDeviation(data)); Assert.Throws(() => Statistics.PopulationVariance(data)); @@ -71,22 +72,27 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests Assert.Throws(() => SortedArrayStatistics.Minimum(data)); Assert.Throws(() => SortedArrayStatistics.Maximum(data)); + Assert.Throws(() => SortedArrayStatistics.OrderStatistic(data, 1)); Assert.Throws(() => SortedArrayStatistics.Median(data)); Assert.Throws(() => SortedArrayStatistics.LowerQuartile(data)); Assert.Throws(() => SortedArrayStatistics.UpperQuartile(data)); Assert.Throws(() => SortedArrayStatistics.Percentile(data, 30)); Assert.Throws(() => SortedArrayStatistics.Quantile(data, 0.3)); - Assert.Throws(() => SortedArrayStatistics.QuantileCompatible(data, 0.3, QuantileCompatibility.Nearest)); + Assert.Throws(() => SortedArrayStatistics.QuantileCustom(data, 0.3, 0, 0, 1, 0)); + Assert.Throws(() => SortedArrayStatistics.QuantileCustom(data, 0.3, QuantileDefinition.Nearest)); Assert.Throws(() => SortedArrayStatistics.InterquartileRange(data)); Assert.Throws(() => SortedArrayStatistics.FiveNumberSummary(data)); Assert.Throws(() => ArrayStatistics.Minimum(data)); Assert.Throws(() => ArrayStatistics.Maximum(data)); + Assert.Throws(() => ArrayStatistics.OrderStatisticInplace(data, 1)); Assert.Throws(() => ArrayStatistics.Mean(data)); Assert.Throws(() => ArrayStatistics.Variance(data)); Assert.Throws(() => ArrayStatistics.StandardDeviation(data)); Assert.Throws(() => ArrayStatistics.PopulationVariance(data)); Assert.Throws(() => ArrayStatistics.PopulationStandardDeviation(data)); + Assert.Throws(() => ArrayStatistics.MedianInplace(data)); + Assert.Throws(() => ArrayStatistics.QuantileInplace(data, 0.3)); Assert.Throws(() => StreamingStatistics.Minimum(data)); Assert.Throws(() => StreamingStatistics.Maximum(data)); @@ -106,6 +112,7 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests Assert.DoesNotThrow(() => Statistics.Maximum(data)); Assert.DoesNotThrow(() => Statistics.Mean(data)); Assert.DoesNotThrow(() => Statistics.Median(data)); + Assert.DoesNotThrow(() => Statistics.Quantile(data, 0.3)); Assert.DoesNotThrow(() => Statistics.Variance(data)); Assert.DoesNotThrow(() => Statistics.StandardDeviation(data)); Assert.DoesNotThrow(() => Statistics.PopulationVariance(data)); @@ -113,22 +120,27 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests Assert.DoesNotThrow(() => SortedArrayStatistics.Minimum(data)); Assert.DoesNotThrow(() => SortedArrayStatistics.Maximum(data)); + Assert.DoesNotThrow(() => SortedArrayStatistics.OrderStatistic(data, 1)); Assert.DoesNotThrow(() => SortedArrayStatistics.Median(data)); Assert.DoesNotThrow(() => SortedArrayStatistics.LowerQuartile(data)); Assert.DoesNotThrow(() => SortedArrayStatistics.UpperQuartile(data)); Assert.DoesNotThrow(() => SortedArrayStatistics.Percentile(data, 30)); Assert.DoesNotThrow(() => SortedArrayStatistics.Quantile(data, 0.3)); - Assert.DoesNotThrow(() => SortedArrayStatistics.QuantileCompatible(data, 0.3, QuantileCompatibility.Nearest)); + Assert.DoesNotThrow(() => SortedArrayStatistics.QuantileCustom(data, 0.3, 0, 0, 1, 0)); + Assert.DoesNotThrow(() => SortedArrayStatistics.QuantileCustom(data, 0.3, QuantileDefinition.Nearest)); Assert.DoesNotThrow(() => SortedArrayStatistics.InterquartileRange(data)); Assert.DoesNotThrow(() => SortedArrayStatistics.FiveNumberSummary(data)); Assert.DoesNotThrow(() => ArrayStatistics.Minimum(data)); Assert.DoesNotThrow(() => ArrayStatistics.Maximum(data)); + Assert.DoesNotThrow(() => ArrayStatistics.OrderStatisticInplace(data, 1)); Assert.DoesNotThrow(() => ArrayStatistics.Mean(data)); Assert.DoesNotThrow(() => ArrayStatistics.Variance(data)); Assert.DoesNotThrow(() => ArrayStatistics.StandardDeviation(data)); Assert.DoesNotThrow(() => ArrayStatistics.PopulationVariance(data)); Assert.DoesNotThrow(() => ArrayStatistics.PopulationStandardDeviation(data)); + Assert.DoesNotThrow(() => ArrayStatistics.MedianInplace(data)); + Assert.DoesNotThrow(() => ArrayStatistics.QuantileInplace(data, 0.3)); Assert.DoesNotThrow(() => StreamingStatistics.Minimum(data)); Assert.DoesNotThrow(() => StreamingStatistics.Maximum(data)); @@ -216,18 +228,292 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests } [Test] - public void MedianOrderOnShortSequence() + public void OrderStatisticsOnShortSequence() { // -3 -1 -0.5 0 1 4 5 6 10 - var samples = new[] {-1, 5, 0, -3, 10, -0.5, 4, 1, 6}; - Assert.That(Statistics.Median(samples), Is.EqualTo(1), "Median"); + var samples = new[] { -1, 5, 0, -3, 10, -0.5, 4, 1, 6 }; + + Assert.That(Statistics.OrderStatistic(samples, 0), Is.NaN, "Order-0 (bad)"); Assert.That(Statistics.OrderStatistic(samples, 1), Is.EqualTo(-3), "Order-1"); + Assert.That(Statistics.OrderStatistic(samples, 2), Is.EqualTo(-1), "Order-2"); Assert.That(Statistics.OrderStatistic(samples, 3), Is.EqualTo(-0.5), "Order-3"); Assert.That(Statistics.OrderStatistic(samples, 7), Is.EqualTo(5), "Order-7"); + Assert.That(Statistics.OrderStatistic(samples, 8), Is.EqualTo(6), "Order-8"); Assert.That(Statistics.OrderStatistic(samples, 9), Is.EqualTo(10), "Order-9"); + Assert.That(Statistics.OrderStatistic(samples, 10), Is.NaN, "Order-10 (bad)"); + + Assert.That(ArrayStatistics.OrderStatisticInplace(samples, 0), Is.NaN, "Order-0 (bad)"); + Assert.That(ArrayStatistics.OrderStatisticInplace(samples, 1), Is.EqualTo(-3), "Order-1"); + Assert.That(ArrayStatistics.OrderStatisticInplace(samples, 2), Is.EqualTo(-1), "Order-2"); + Assert.That(ArrayStatistics.OrderStatisticInplace(samples, 3), Is.EqualTo(-0.5), "Order-3"); + Assert.That(ArrayStatistics.OrderStatisticInplace(samples, 7), Is.EqualTo(5), "Order-7"); + Assert.That(ArrayStatistics.OrderStatisticInplace(samples, 8), Is.EqualTo(6), "Order-8"); + Assert.That(ArrayStatistics.OrderStatisticInplace(samples, 9), Is.EqualTo(10), "Order-9"); + Assert.That(ArrayStatistics.OrderStatisticInplace(samples, 10), Is.NaN, "Order-10 (bad)"); + + Array.Sort(samples); + Assert.That(SortedArrayStatistics.OrderStatistic(samples, 0), Is.NaN, "Order-0 (bad)"); + Assert.That(SortedArrayStatistics.OrderStatistic(samples, 1), Is.EqualTo(-3), "Order-1"); + Assert.That(SortedArrayStatistics.OrderStatistic(samples, 2), Is.EqualTo(-1), "Order-2"); + Assert.That(SortedArrayStatistics.OrderStatistic(samples, 3), Is.EqualTo(-0.5), "Order-3"); + Assert.That(SortedArrayStatistics.OrderStatistic(samples, 7), Is.EqualTo(5), "Order-7"); + Assert.That(SortedArrayStatistics.OrderStatistic(samples, 8), Is.EqualTo(6), "Order-8"); + Assert.That(SortedArrayStatistics.OrderStatistic(samples, 9), Is.EqualTo(10), "Order-9"); + Assert.That(SortedArrayStatistics.OrderStatistic(samples, 10), Is.NaN, "Order-10 (bad)"); + } + + [TestCase(0d, -3d)] + [TestCase(1d, 10d)] + [TestCase(0.5d, 1/5d)] + [TestCase(0.2d, -1d)] + [TestCase(0.7d, 4d)] + [TestCase(0.01d, -3d)] + [TestCase(0.99d, 10d)] + [TestCase(0.52d, 1d)] + [TestCase(0.325d, 0d)] + public void QuantileR1InverseCDFOnShortSequence(double tau, double expected) + { + // R: quantile(c(-1,5,0,-3,10,-0.5,4,0.2,1,6),probs=c(0,1,0.5,0.2,0.7,0.01,0.99,0.52,0.325),type=1) + // Mathematica: Quantile[{-1,5,0,-3,10,-1/2,4,1/5,1,6},{0,1,1/2,1/5,7/10,1/100,99/100,13/25,13/40},{{0,0},{1,0}}] + + var samples = new[] {-1, 5, 0, -3, 10, -0.5, 4, 0.2, 1, 6}; + + Assert.AreEqual(expected, Statistics.InverseCDF(samples, tau), 1e-14); + Assert.AreEqual(expected, Statistics.QuantileCustom(samples, tau, QuantileDefinition.InverseCDF), 1e-14); + + Assert.AreEqual(expected, ArrayStatistics.QuantileCustomInplace(samples, tau, QuantileDefinition.InverseCDF), 1e-14); + Assert.AreEqual(expected, ArrayStatistics.QuantileCustomInplace(samples, tau, 0d, 0d, 1d, 0d), 1e-14); + + Array.Sort(samples); + Assert.AreEqual(expected, SortedArrayStatistics.QuantileCustom(samples, tau, QuantileDefinition.InverseCDF), 1e-14); + Assert.AreEqual(expected, SortedArrayStatistics.QuantileCustom(samples, tau, 0d, 0d, 1d, 0d), 1e-14); + } + + [TestCase(0d, -3d)] + [TestCase(1d, 10d)] + [TestCase(0.5d, 3/5d)] + [TestCase(0.2d, -3/4d)] + [TestCase(0.7d, 9/2d)] + [TestCase(0.01d, -3d)] + [TestCase(0.99d, 10d)] + [TestCase(0.52d, 1d)] + [TestCase(0.325d, 0d)] + public void QuantileR2InverseCDFAverageOnShortSequence(double tau, double expected) + { + // R: quantile(c(-1,5,0,-3,10,-0.5,4,0.2,1,6),probs=c(0,1,0.5,0.2,0.7,0.01,0.99,0.52,0.325),type=2) + // Mathematica: Not Supported + + var samples = new[] {-1, 5, 0, -3, 10, -0.5, 4, 0.2, 1, 6}; + + Assert.AreEqual(expected, Statistics.QuantileCustom(samples, tau, QuantileDefinition.R2), 1e-14); + + Assert.AreEqual(expected, ArrayStatistics.QuantileCustomInplace(samples, tau, QuantileDefinition.InverseCDFAverage), 1e-14); + + Array.Sort(samples); + Assert.AreEqual(expected, SortedArrayStatistics.QuantileCustom(samples, tau, QuantileDefinition.InverseCDFAverage), 1e-14); + } + + [TestCase(0d, -3d)] + [TestCase(1d, 10d)] + [TestCase(0.5d, 1/5d)] + [TestCase(0.2d, -1d)] + [TestCase(0.7d, 4d)] + [TestCase(0.01d, -3d)] + [TestCase(0.99d, 10d)] + [TestCase(0.52d, 1/5d)] + [TestCase(0.325d, -1/2d)] + public void QuantileR3NearestOnShortSequence(double tau, double expected) + { + // R: quantile(c(-1,5,0,-3,10,-0.5,4,0.2,1,6),probs=c(0,1,0.5,0.2,0.7,0.01,0.99,0.52,0.325),type=3) + // Mathematica: Quantile[{-1,5,0,-3,10,-1/2,4,1/5,1,6},{0,1,1/2,1/5,7/10,1/100,99/100,13/25,13/40},{{1/2,0},{0,0}}] + + var samples = new[] {-1, 5, 0, -3, 10, -0.5, 4, 0.2, 1, 6}; + + Assert.AreEqual(expected, Statistics.QuantileCustom(samples, tau, QuantileDefinition.R3), 1e-14); + + Assert.AreEqual(expected, ArrayStatistics.QuantileCustomInplace(samples, tau, QuantileDefinition.Nearest), 1e-14); + Assert.AreEqual(expected, ArrayStatistics.QuantileCustomInplace(samples, tau, 0.5d, 0d, 0d, 0d), 1e-14); + + Array.Sort(samples); + Assert.AreEqual(expected, SortedArrayStatistics.QuantileCustom(samples, tau, QuantileDefinition.Nearest), 1e-14); + Assert.AreEqual(expected, SortedArrayStatistics.QuantileCustom(samples, tau, 0.5d, 0d, 0d, 0d), 1e-14); + } + + [TestCase(0d, -3d)] + [TestCase(1d, 10d)] + [TestCase(0.5d, 1/5d)] + [TestCase(0.2d, -1d)] + [TestCase(0.7d, 4d)] + [TestCase(0.01d, -3d)] + [TestCase(0.99d, 48/5d)] + [TestCase(0.52d, 9/25d)] + [TestCase(0.325d, -3/8d)] + public void QuantileR4CaliforniaOnShortSequence(double tau, double expected) + { + // R: quantile(c(-1,5,0,-3,10,-0.5,4,0.2,1,6),probs=c(0,1,0.5,0.2,0.7,0.01,0.99,0.52,0.325),type=4) + // Mathematica: Quantile[{-1,5,0,-3,10,-1/2,4,1/5,1,6},{0,1,1/2,1/5,7/10,1/100,99/100,13/25,13/40},{{0,0},{0,1}}] + + var samples = new[] {-1, 5, 0, -3, 10, -0.5, 4, 0.2, 1, 6}; + + Assert.AreEqual(expected, Statistics.QuantileCustom(samples, tau, QuantileDefinition.R4), 1e-14); + + Assert.AreEqual(expected, ArrayStatistics.QuantileCustomInplace(samples, tau, QuantileDefinition.California), 1e-14); + Assert.AreEqual(expected, ArrayStatistics.QuantileCustomInplace(samples, tau, 0d, 0d, 0d, 1d), 1e-14); + + Array.Sort(samples); + Assert.AreEqual(expected, SortedArrayStatistics.QuantileCustom(samples, tau, QuantileDefinition.California), 1e-14); + Assert.AreEqual(expected, SortedArrayStatistics.QuantileCustom(samples, tau, 0d, 0d, 0d, 1d), 1e-14); + } + + [TestCase(0d, -3d)] + [TestCase(1d, 10d)] + [TestCase(0.5d, 3/5d)] + [TestCase(0.2d, -3/4d)] + [TestCase(0.7d, 9/2d)] + [TestCase(0.01d, -3d)] + [TestCase(0.99d, 10d)] + [TestCase(0.52d, 19/25d)] + [TestCase(0.325d, -1/8d)] + public void QuantileR5HydrologyOnShortSequence(double tau, double expected) + { + // R: quantile(c(-1,5,0,-3,10,-0.5,4,0.2,1,6),probs=c(0,1,0.5,0.2,0.7,0.01,0.99,0.52,0.325),type=5) + // Mathematica: Quantile[{-1,5,0,-3,10,-1/2,4,1/5,1,6},{0,1,1/2,1/5,7/10,1/100,99/100,13/25,13/40},{{1/2,0},{0,1}}] + + var samples = new[] {-1, 5, 0, -3, 10, -0.5, 4, 0.2, 1, 6}; + + Assert.AreEqual(expected, Statistics.QuantileCustom(samples, tau, QuantileDefinition.R5), 1e-14); + + Assert.AreEqual(expected, ArrayStatistics.QuantileCustomInplace(samples, tau, QuantileDefinition.Hydrology), 1e-14); + Assert.AreEqual(expected, ArrayStatistics.QuantileCustomInplace(samples, tau, 0.5d, 0d, 0d, 1d), 1e-14); + + Array.Sort(samples); + Assert.AreEqual(expected, SortedArrayStatistics.QuantileCustom(samples, tau, QuantileDefinition.Hydrology), 1e-14); + Assert.AreEqual(expected, SortedArrayStatistics.QuantileCustom(samples, tau, 0.5d, 0d, 0d, 1d), 1e-14); + } + + [TestCase(0d, -3d)] + [TestCase(1d, 10d)] + [TestCase(0.5d, 3/5d)] + [TestCase(0.2d, -9/10d)] + [TestCase(0.7d, 47/10d)] + [TestCase(0.01d, -3d)] + [TestCase(0.99d, 10d)] + [TestCase(0.52d, 97/125d)] + [TestCase(0.325d, -17/80d)] + public void QuantileR6WeibullOnShortSequence(double tau, double expected) + { + // R: quantile(c(-1,5,0,-3,10,-0.5,4,0.2,1,6),probs=c(0,1,0.5,0.2,0.7,0.01,0.99,0.52,0.325),type=6) + // Mathematica: Quantile[{-1,5,0,-3,10,-1/2,4,1/5,1,6},{0,1,1/2,1/5,7/10,1/100,99/100,13/25,13/40},{{0,1},{0,1}}] + + var samples = new[] {-1, 5, 0, -3, 10, -0.5, 4, 0.2, 1, 6}; + + Assert.AreEqual(expected, Statistics.QuantileCustom(samples, tau, QuantileDefinition.R6), 1e-14); + + Assert.AreEqual(expected, ArrayStatistics.QuantileCustomInplace(samples, tau, QuantileDefinition.Weibull), 1e-14); + Assert.AreEqual(expected, ArrayStatistics.QuantileCustomInplace(samples, tau, 0d, 1d, 0d, 1d), 1e-14); + + Array.Sort(samples); + Assert.AreEqual(expected, SortedArrayStatistics.QuantileCustom(samples, tau, QuantileDefinition.Weibull), 1e-14); + Assert.AreEqual(expected, SortedArrayStatistics.QuantileCustom(samples, tau, 0d, 1d, 0d, 1d), 1e-14); + } + + [TestCase(0d, -3d)] + [TestCase(1d, 10d)] + [TestCase(0.5d, 3/5d)] + [TestCase(0.2d, -3/5d)] + [TestCase(0.7d, 43/10d)] + [TestCase(0.01d, -141/50d)] + [TestCase(0.99d, 241/25d)] + [TestCase(0.52d, 93/125d)] + [TestCase(0.325d, -3/80d)] + public void QuantileR7ExcelOnShortSequence(double tau, double expected) + { + // R: quantile(c(-1,5,0,-3,10,-0.5,4,0.2,1,6),probs=c(0,1,0.5,0.2,0.7,0.01,0.99,0.52,0.325),type=7) + // Mathematica: Quantile[{-1,5,0,-3,10,-1/2,4,1/5,1,6},{0,1,1/2,1/5,7/10,1/100,99/100,13/25,13/40},{{1,-1},{0,1}}] + + var samples = new[] {-1, 5, 0, -3, 10, -0.5, 4, 0.2, 1, 6}; + + Assert.AreEqual(expected, Statistics.QuantileCustom(samples, tau, QuantileDefinition.R7), 1e-14); + + Assert.AreEqual(expected, ArrayStatistics.QuantileCustomInplace(samples, tau, QuantileDefinition.Excel), 1e-14); + Assert.AreEqual(expected, ArrayStatistics.QuantileCustomInplace(samples, tau, 1d, -1d, 0d, 1d), 1e-14); Array.Sort(samples); - Assert.That(SortedArrayStatistics.Median(samples), Is.EqualTo(1), "Median"); + Assert.AreEqual(expected, SortedArrayStatistics.QuantileCustom(samples, tau, QuantileDefinition.Excel), 1e-14); + Assert.AreEqual(expected, SortedArrayStatistics.QuantileCustom(samples, tau, 1d, -1d, 0d, 1d), 1e-14); + } + + [TestCase(0d, -3d)] + [TestCase(1d, 10d)] + [TestCase(0.5d, 3/5d)] + [TestCase(0.2d, -4/5d)] + [TestCase(0.7d, 137/30d)] + [TestCase(0.01d, -3d)] + [TestCase(0.99d, 10d)] + [TestCase(0.52d, 287/375d)] + [TestCase(0.325d, -37/240d)] + public void QuantileR8MedianOnShortSequence(double tau, double expected) + { + // R: quantile(c(-1,5,0,-3,10,-0.5,4,0.2,1,6),probs=c(0,1,0.5,0.2,0.7,0.01,0.99,0.52,0.325),type=8) + // Mathematica: Quantile[{-1,5,0,-3,10,-1/2,4,1/5,1,6},{0,1,1/2,1/5,7/10,1/100,99/100,13/25,13/40},{{1/3,1/3},{0,1}}] + var samples = new[] {-1, 5, 0, -3, 10, -0.5, 4, 0.2, 1, 6}; + + Assert.AreEqual(expected, Statistics.Quantile(samples, tau), 1e-14); + Assert.AreEqual(expected, Statistics.QuantileCustom(samples, tau, QuantileDefinition.R8), 1e-14); + + Assert.AreEqual(expected, ArrayStatistics.QuantileInplace(samples, tau), 1e-14); + Assert.AreEqual(expected, ArrayStatistics.QuantileCustomInplace(samples, tau, QuantileDefinition.Median), 1e-14); + Assert.AreEqual(expected, ArrayStatistics.QuantileCustomInplace(samples, tau, 1 / 3d, 1 / 3d, 0d, 1d), 1e-14); + + Array.Sort(samples); + Assert.AreEqual(expected, SortedArrayStatistics.Quantile(samples, tau), 1e-14); + Assert.AreEqual(expected, SortedArrayStatistics.QuantileCustom(samples, tau, QuantileDefinition.Median), 1e-14); + Assert.AreEqual(expected, SortedArrayStatistics.QuantileCustom(samples, tau, 1/3d, 1/3d, 0d, 1d), 1e-14); + } + + [TestCase(0d, -3d)] + [TestCase(1d, 10d)] + [TestCase(0.5d, 3/5d)] + [TestCase(0.2d, -63/80d)] + [TestCase(0.7d, 91/20d)] + [TestCase(0.01d, -3d)] + [TestCase(0.99d, 10d)] + [TestCase(0.52d, 191/250d)] + [TestCase(0.325d, -47/320d)] + public void QuantileR9NormalOnShortSequence(double tau, double expected) + { + // R: quantile(c(-1,5,0,-3,10,-0.5,4,0.2,1,6),probs=c(0,1,0.5,0.2,0.7,0.01,0.99,0.52,0.325),type=9) + // Mathematica: Quantile[{-1,5,0,-3,10,-1/2,4,1/5,1,6},{0,1,1/2,1/5,7/10,1/100,99/100,13/25,13/40},{{3/8,1/4},{0,1}}] + var samples = new[] {-1, 5, 0, -3, 10, -0.5, 4, 0.2, 1, 6}; + + Assert.AreEqual(expected, Statistics.QuantileCustom(samples, tau, QuantileDefinition.R9), 1e-14); + + Assert.AreEqual(expected, ArrayStatistics.QuantileCustomInplace(samples, tau, QuantileDefinition.Normal), 1e-14); + Assert.AreEqual(expected, ArrayStatistics.QuantileCustomInplace(samples, tau, 3/8d, 1/4d, 0d, 1d), 1e-14); + + Array.Sort(samples); + Assert.AreEqual(expected, SortedArrayStatistics.QuantileCustom(samples, tau, QuantileDefinition.Normal), 1e-14); + Assert.AreEqual(expected, SortedArrayStatistics.QuantileCustom(samples, tau, 3/8d, 1/4d, 0d, 1d), 1e-14); + } + + [Test] + public void MedianOnShortSequence() + { + // R: median(c(-1,5,0,-3,10,-0.5,4,0.2,1,6)) + // Mathematica: Median[{-1,5,0,-3,10,-1/2,4,1/5,1,6}] + var even = new[] {-1, 5, 0, -3, 10, -0.5, 4, 0.2, 1, 6}; + Assert.AreEqual(0.6d, Statistics.Median(even), 1e-14); + Assert.AreEqual(0.6d, ArrayStatistics.MedianInplace(even), 1e-14); + Array.Sort(even); + Assert.AreEqual(0.6d, SortedArrayStatistics.Median(even), 1e-14); + + // R: median(c(-1,5,0,-3,10,-0.5,4,0.2,1)) + // Mathematica: Median[{-1,5,0,-3,10,-1/2,4,1/5,1}] + var odd = new[] { -1, 5, 0, -3, 10, -0.5, 4, 0.2, 1 }; + Assert.AreEqual(0.2d, Statistics.Median(odd), 1e-14); + Assert.AreEqual(0.2d, ArrayStatistics.MedianInplace(odd), 1e-14); + Array.Sort(even); + Assert.AreEqual(0.2d, SortedArrayStatistics.Median(odd), 1e-14); } /// @@ -331,9 +617,11 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests var seq = File.ReadLines("./data/Codeplex-5667.csv").Select(double.Parse); Assert.AreEqual(1.0, Statistics.Median(seq)); - var sorted = seq.ToArray(); - Array.Sort(sorted); - Assert.AreEqual(1.0, SortedArrayStatistics.Median(sorted)); + var array = seq.ToArray(); + Assert.AreEqual(1.0, ArrayStatistics.MedianInplace(array)); + + Array.Sort(array); + Assert.AreEqual(1.0, SortedArrayStatistics.Median(array)); } } }