Statistics: Proper Quantile statistics (supporting all 9 R-types and custom)

src/Numerics/Numerics.csproj

@@ -112,6 +112,7 @@
     <Compile Include="SpecialFunctions\ModifiedBessel.cs" />
     <Compile Include="SpecialFunctions\Logistic.cs" />
     <Compile Include="Statistics\ArrayStatistics.cs" />
+    <Compile Include="Statistics\QuantileDefinition.cs" />
     <Compile Include="Statistics\StreamingStatistics.cs" />
     <Compile Include="Statistics\SortedArrayStatistics.cs" />
     <Compile Include="TargetedPatchingOptOutAttribute.cs" />

src/Numerics/Statistics/ArrayStatistics.cs

@@ -32,6 +32,13 @@ using System;
 
 namespace MathNet.Numerics.Statistics
 {
+    /// <summary>
+    /// Statistics operating on arrays assumed to be unsorted.
+    /// WARNING: Methods with the Inplace-suffix may modify the data array by reordering its entries.
+    /// </summary>
+    /// <seealso cref="SortedArrayStatistics"/>
+    /// <seealso cref="StreamingStatistics"/>
+    /// <seealso cref="Statistics"/>
     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;
         }
 
+        /// <summary>
+        /// 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.
+        /// </summary>
+        /// <param name="data">Sample array, no sorting is assumed. Will be reordered.</param>
+        /// <param name="order">One-based order of the statistic, must be between 1 and N (inclusive).</param>
+        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);
+        }
+
         /// <summary>
         /// 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));
         }
+
+        /// <summary>
+        /// 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.
+        /// </summary>
+        /// <param name="data">Sample array, no sorting is assumed. Will be reordered.</param>
+        public static double MedianInplace(double[] data)
+        {
+            return QuantileInplace(data, 0.5d);
+        }
+
+        /// <summary>
+        /// 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.
+        /// </summary>
+        /// <param name="data">Sample array, no sorting is assumed. Will be reordered.</param>
+        /// <param name="tau">Quantile selector, between 0.0 and 1.0 (inclusive).</param>
+        /// <remarks>
+        /// R-8, SciPy-(1/3,1/3):
+        /// Linear interpolation of the approximate medians for order statistics.
+        /// When tau &lt; (2/3) / (N + 1/3), use x1. When tau &gt;= (N - 1/3) / (N + 1/3), use xN.
+        /// </remarks>
+        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);
+        }
+
+        /// <summary>
+        /// 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.
+        /// </summary>
+        /// <param name="data">Sample array, no sorting is assumed. Will be reordered.</param>
+        /// <param name="tau">Quantile selector, between 0.0 and 1.0 (inclusive)</param>
+        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);
+        }
+
+        /// <summary>
+        /// 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.
+        /// </summary>
+        /// <param name="data">Sample array, no sorting is assumed. Will be reordered.</param>
+        /// <param name="tau">Quantile selector, between 0.0 and 1.0 (inclusive)</param>
+        /// <param name="definition">Quantile definition, to choose what product/definition it should be consistent with</param>
+        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;
+            }
+        }
     }
 }

src/Numerics/Statistics/Histogram.cs

@@ -246,7 +246,7 @@ namespace MathNet.Numerics.Statistics
         /// </summary>
         /// <param name="data">The datasequence to build a histogram on.</param>
         /// <param name="nbuckets">The number of buckets to use.</param>
-        public Histogram(IEnumerable<double> data, int nbuckets) 
+        public Histogram(IEnumerable<double> 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.

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);
             }

src/Numerics/Statistics/QuantileDefinition.cs

@@ -0,0 +1,45 @@
+// <copyright file="QuantileDefinition.cs" company="Math.NET">
+// 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.
+// </copyright>
+
+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,
+    }
+}

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
-    }
-
+    /// <summary>
+    /// Statistics operating on an array already sorted ascendingly.
+    /// </summary>
+    /// <seealso cref="ArrayStatistics"/>
+    /// <seealso cref="StreamingStatistics"/>
+    /// <seealso cref="Statistics"/>
     public static class SortedArrayStatistics
     {
-        const double Third = 1d/3d;
-        const double Half = 1d/2d;
-
         /// <summary>
         /// Returns the smallest value from the sorted data array (ascending).
         /// </summary>
@@ -69,6 +64,19 @@ namespace MathNet.Numerics.Statistics
             return data[data.Length - 1];
         }
 
+        /// <summary>
+        /// Returns the order statistic (order 1..N) from the sorted data array (ascending).
+        /// </summary>
+        /// <param name="data">Sample array, must be sorted ascendingly.</param>
+        /// <param name="order">One-based order of the statistic, must be between 1 and N (inclusive).</param>
+        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];
+        }
+
         /// <summary>
         /// 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]);
+        }
+
+        /// <summary>
+        /// 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.
+        /// </summary>
+        /// <param name="data">Sample array, must be sorted ascendingly.</param>
+        /// <param name="tau">Quantile selector, between 0.0 and 1.0 (inclusive).</param>
+        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);
         }
 
         /// <summary>
         /// 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.
         /// </summary>
-        public static double QuantileCompatible(double[] data, double tau, QuantileCompatibility compatibility)
+        /// <param name="data">Sample array, must be sorted ascendingly.</param>
+        /// <param name="tau">Quantile selector, between 0.0 and 1.0 (inclusive).</param>
+        /// <param name="definition">Quantile definition, to choose what product/definition it should be consistent with</param>
+        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();

src/Numerics/Statistics/Statistics.cs

@@ -33,7 +33,6 @@ namespace MathNet.Numerics.Statistics
     using System;
     using System.Collections.Generic;
     using System.Linq;
-    using Properties;
 
     /// <summary>
     /// Extension methods to return basic statistics on set of data.
@@ -88,7 +87,7 @@ namespace MathNet.Numerics.Statistics
         }
 
         /// <summary>
-        /// Calculates the sample mean.
+        /// Estimates the sample mean.
         /// </summary>
         /// <param name="data">The data to calculate the mean of.</param>
         /// <returns>The mean of the sample.</returns>
@@ -101,7 +100,7 @@ namespace MathNet.Numerics.Statistics
         }
 
         /// <summary>
-        /// Calculates the sample mean.
+        /// Estimates the sample mean.
         /// </summary>
         /// <param name="data">The data to calculate the mean of.</param>
         /// <returns>The mean of the sample.</returns>
@@ -112,7 +111,7 @@ namespace MathNet.Numerics.Statistics
         }
 
         /// <summary>
-        /// 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).
         /// </summary>
         /// <param name="data">The data to calculate the variance of.</param>
         /// <returns>The unbiased population variance of the sample.</returns>
@@ -125,7 +124,7 @@ namespace MathNet.Numerics.Statistics
         }
 
         /// <summary>
-        /// 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.
         /// </summary>
         /// <param name="data">The data to calculate the variance of.</param>
         /// <returns>The population variance of the sample.</returns>
@@ -136,7 +135,7 @@ namespace MathNet.Numerics.Statistics
         }
 
         /// <summary>
-        /// 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).
         /// </summary>
         /// <param name="data">The data to calculate the variance of.</param>
         /// <returns>The biased population variance of the sample.</returns>
@@ -149,7 +148,7 @@ namespace MathNet.Numerics.Statistics
         }
 
         /// <summary>
-        /// 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.
         /// </summary>
         /// <param name="data">The data to calculate the variance of.</param>
         /// <returns>The population variance of the sample.</returns>
@@ -160,7 +159,7 @@ namespace MathNet.Numerics.Statistics
         }
 
         /// <summary>
-        /// 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).
         /// </summary>
         /// <param name="data">The data to calculate the standard deviation of.</param>
         /// <returns>The standard deviation of the sample.</returns>
@@ -173,7 +172,7 @@ namespace MathNet.Numerics.Statistics
         }
 
         /// <summary>
-        /// 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).
         /// </summary>
         /// <param name="data">The data to calculate the standard deviation of.</param>
         /// <returns>The standard deviation of the sample.</returns>
@@ -184,7 +183,7 @@ namespace MathNet.Numerics.Statistics
         }
 
         /// <summary>
-        /// 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).
         /// </summary>
         /// <param name="data">The data to calculate the standard deviation of.</param>
         /// <returns>The standard deviation of the sample.</returns>
@@ -197,7 +196,7 @@ namespace MathNet.Numerics.Statistics
         }
 
         /// <summary>
-        /// 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).
         /// </summary>
         /// <param name="data">The data to calculate the standard deviation of.</param>
         /// <returns>The standard deviation of the sample.</returns>
@@ -208,158 +207,120 @@ namespace MathNet.Numerics.Statistics
         }
 
         /// <summary>
-        /// Calculates the sample median.
+        /// Estimates the sample median.
         /// </summary>
         /// <param name="data">The data to calculate the median of.</param>
         /// <returns>The median of the sample.</returns>
         public static double Median(this IEnumerable<double> data)
         {
-            if (data == null)
-            {
-                throw new ArgumentNullException("data");
-            }
-
-            var dataArray = new List<double>(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);
         }
 
         /// <summary>
-        /// Calculates the sample median.
+        /// Estimates the sample median.
         /// </summary>
         /// <param name="data">The data to calculate the median of.</param>
         /// <returns>The median of the sample.</returns>
         public static double Median(this IEnumerable<double?> data)
         {
-            if (data == null)
-            {
-                throw new ArgumentNullException("data");
-            }
-
-            var nonNull = new List<double>();
-            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);
         }
 
         /// <summary>
-        /// Evaluate the i-order (1..N) statistic of the provided samples.
+        /// Estimates the sample tau-quantile.
         /// </summary>
-        /// <param name="samples">The sample data.</param>
-        /// <param name="order">Order of the statistic to evaluate.</param>
-        /// <returns>The i'th order statistic in the sample data.</returns>
-        public static double OrderStatistic(IEnumerable<double> samples, int order)
+        /// <param name="data">The data to calculate the median of.</param>
+        /// <param name="tau">Quantile selector, between 0.0 and 1.0 (inclusive).</param>
+        /// <returns>The median of the sample.</returns>
+        public static double Quantile(this IEnumerable<double> data, double tau)
         {
-            if (order == 1)
-            {
-                // Can be done in linear time by Min()
-                return Minimum(samples);
-            }
-
-            var list = new List<double>(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);
         }
 
         /// <summary>
-        /// 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.
         /// </summary>
-        /// <param name="samples">The sample data.</param>
-        /// <param name="left">The left bound in which to order select.</param>
-        /// <param name="right">The right bound in which to order select.</param>
-        /// <param name="order">The order we are trying to find.</param>
-        /// <returns>The <paramref name="order"/> order statistic.</returns>
-        private static double OrderSelect(IList<double> samples, int left, int right, int order)
+        /// <param name="data">The data to calculate the median of.</param>
+        /// <param name="tau">Quantile selector, between 0.0 and 1.0 (inclusive).</param>
+        /// <returns>The median of the sample.</returns>
+        public static double Quantile(this IEnumerable<double?> 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];
+        /// <summary>
+        /// Estimates the empiric inverse CDF at tau (tau-quantile).
+        /// </summary>
+        /// <param name="data">The data to calculate the median of.</param>
+        /// <param name="tau">Quantile selector, between 0.0 and 1.0 (inclusive).</param>
+        /// <returns>The median of the sample.</returns>
+        public static double InverseCDF(this IEnumerable<double> 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);
-                    }
-                }
+        /// <summary>
+        /// Estimates the empiric inverse CDF at tau (tau-quantile).
+        /// </summary>
+        /// <param name="data">The data to calculate the median of.</param>
+        /// <param name="tau">Quantile selector, between 0.0 and 1.0 (inclusive).</param>
+        /// <returns>The median of the sample.</returns>
+        public static double InverseCDF(this IEnumerable<double?> 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);
+        /// <summary>
+        /// Estimates the sample tau-quantile.
+        /// </summary>
+        /// <param name="data">The data to calculate the median of.</param>
+        /// <param name="tau">Quantile selector, between 0.0 and 1.0 (inclusive).</param>
+        /// <returns>The median of the sample.</returns>
+        /// <param name="definition">Quantile definition, to choose what product/definition it should be consistent with</param>
+        public static double QuantileCustom(this IEnumerable<double> 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;
-                }
+        /// <summary>
+        /// Estimates the sample tau-quantile.
+        /// </summary>
+        /// <param name="data">The data to calculate the median of.</param>
+        /// <param name="tau">Quantile selector, between 0.0 and 1.0 (inclusive).</param>
+        /// <returns>The median of the sample.</returns>
+        /// <param name="definition">Quantile definition, to choose what product/definition it should be consistent with</param>
+        public static double QuantileCustom(this IEnumerable<double?> 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;
-                }
-            }
+        /// <summary>
+        /// Returns the i-order (1..N) statistic of the provided samples.
+        /// </summary>
+        /// <param name="data">The sample data.</param>
+        /// <param name="order">Order of the statistic to evaluate.</param>
+        /// <returns>The i'th order statistic in the sample data.</returns>
+        public static double OrderStatistic(IEnumerable<double> data, int order)
+        {
+            if (data == null) throw new ArgumentNullException("data");
+            var array = data.ToArray();
+            return ArrayStatistics.OrderStatisticInplace(array, order);
         }
     }
 }

src/Numerics/Statistics/StreamingStatistics.cs

@@ -33,6 +33,13 @@ using System.Collections.Generic;
 
 namespace MathNet.Numerics.Statistics
 {
+    /// <summary>
+    /// 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.
+    /// </summary>
+    /// <seealso cref="SortedArrayStatistics"/>
+    /// <seealso cref="StreamingStatistics"/>
+    /// <seealso cref="Statistics"/>
     public static class StreamingStatistics
     {
         /// <summary>

src/Portable/Portable.csproj

@@ -1062,6 +1062,9 @@
     <Compile Include="..\Numerics\Statistics\Percentile.cs">
       <Link>Statistics\Percentile.cs</Link>
     </Compile>
+    <Compile Include="..\Numerics\Statistics\QuantileDefinition.cs">
+      <Link>Statistics\QuantileDefinition.cs</Link>
+    </Compile>
     <Compile Include="..\Numerics\Statistics\SortedArrayStatistics.cs">
       <Link>Statistics\SortedArrayStatistics.cs</Link>
     </Compile>

src/UnitTests/StatisticsTests/StatisticsTests.cs

@@ -64,6 +64,7 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests
             Assert.Throws<ArgumentNullException>(() => Statistics.Maximum(data));
             Assert.Throws<ArgumentNullException>(() => Statistics.Mean(data));
             Assert.Throws<ArgumentNullException>(() => Statistics.Median(data));
+            Assert.Throws<ArgumentNullException>(() => Statistics.Quantile(data, 0.3));
             Assert.Throws<ArgumentNullException>(() => Statistics.Variance(data));
             Assert.Throws<ArgumentNullException>(() => Statistics.StandardDeviation(data));
             Assert.Throws<ArgumentNullException>(() => Statistics.PopulationVariance(data));
@@ -71,22 +72,27 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests
 
             Assert.Throws<ArgumentNullException>(() => SortedArrayStatistics.Minimum(data));
             Assert.Throws<ArgumentNullException>(() => SortedArrayStatistics.Maximum(data));
+            Assert.Throws<ArgumentNullException>(() => SortedArrayStatistics.OrderStatistic(data, 1));
             Assert.Throws<ArgumentNullException>(() => SortedArrayStatistics.Median(data));
             Assert.Throws<ArgumentNullException>(() => SortedArrayStatistics.LowerQuartile(data));
             Assert.Throws<ArgumentNullException>(() => SortedArrayStatistics.UpperQuartile(data));
             Assert.Throws<ArgumentNullException>(() => SortedArrayStatistics.Percentile(data, 30));
             Assert.Throws<ArgumentNullException>(() => SortedArrayStatistics.Quantile(data, 0.3));
-            Assert.Throws<ArgumentNullException>(() => SortedArrayStatistics.QuantileCompatible(data, 0.3, QuantileCompatibility.Nearest));
+            Assert.Throws<ArgumentNullException>(() => SortedArrayStatistics.QuantileCustom(data, 0.3, 0, 0, 1, 0));
+            Assert.Throws<ArgumentNullException>(() => SortedArrayStatistics.QuantileCustom(data, 0.3, QuantileDefinition.Nearest));
             Assert.Throws<ArgumentNullException>(() => SortedArrayStatistics.InterquartileRange(data));
             Assert.Throws<ArgumentNullException>(() => SortedArrayStatistics.FiveNumberSummary(data));
 
             Assert.Throws<ArgumentNullException>(() => ArrayStatistics.Minimum(data));
             Assert.Throws<ArgumentNullException>(() => ArrayStatistics.Maximum(data));
+            Assert.Throws<ArgumentNullException>(() => ArrayStatistics.OrderStatisticInplace(data, 1));
             Assert.Throws<ArgumentNullException>(() => ArrayStatistics.Mean(data));
             Assert.Throws<ArgumentNullException>(() => ArrayStatistics.Variance(data));
             Assert.Throws<ArgumentNullException>(() => ArrayStatistics.StandardDeviation(data));
             Assert.Throws<ArgumentNullException>(() => ArrayStatistics.PopulationVariance(data));
             Assert.Throws<ArgumentNullException>(() => ArrayStatistics.PopulationStandardDeviation(data));
+            Assert.Throws<ArgumentNullException>(() => ArrayStatistics.MedianInplace(data));
+            Assert.Throws<ArgumentNullException>(() => ArrayStatistics.QuantileInplace(data, 0.3));
 
             Assert.Throws<ArgumentNullException>(() => StreamingStatistics.Minimum(data));
             Assert.Throws<ArgumentNullException>(() => 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);
         }
 
         /// <summary>
@@ -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));
         }
     }
 }