Statistics: introduce SortedArrayStatistics, use it in Percentile class

src/Numerics/Numerics.csproj

@@ -111,6 +111,7 @@
     <Compile Include="SpecialFunctions\ModifiedStruve.cs" />
     <Compile Include="SpecialFunctions\ModifiedBessel.cs" />
     <Compile Include="SpecialFunctions\Logistic.cs" />
+    <Compile Include="Statistics\SortedArrayStatistics.cs" />
     <Compile Include="TargetedPatchingOptOutAttribute.cs" />
     <Compile Include="Distributions\Continuous\Cauchy.cs" />
     <Compile Include="Distributions\Continuous\Chi.cs" />

src/Numerics/Statistics/Percentile.cs

@@ -71,7 +71,7 @@ namespace MathNet.Numerics.Statistics
         /// <summary>
         /// Holds the data.
         /// </summary>
-        private readonly List<double> _data;
+        private readonly double[] _data;
 
         /// <summary>
         /// Gets or sets the method used to calculate the percentiles.
@@ -94,9 +94,9 @@ namespace MathNet.Numerics.Statistics
             {
                 throw new ArgumentNullException("data");
             }
-            
+
-            _data = new List<double>(data);
+            _data = data.ToArray();
-            _data.Sort();
+            Array.Sort(_data);
         }
 
         /// <summary>
@@ -106,39 +106,19 @@ namespace MathNet.Numerics.Statistics
         /// <returns>the requested percentile.</returns>
         public double Compute(double percentile)
         {
-            if (percentile < 0 || percentile > 1 || _data.Count == 0)
-            {
-                return double.NaN;
-            }
-
-            if (percentile == 0.0 || _data.Count == 1)
-            {
-                return _data[0];
-            }
-
-            if (percentile == 1.0)
-            {
-                return _data[_data.Count - 1];
-            }
-
-            var result = double.NaN;
             switch (Method)
             {
                 case PercentileMethod.Nist:
-                    result = Nist(percentile);
+                    return SortedArrayStatistics.QuantileCompatible(_data, percentile, QuantileCompatibility.Nist);
-                    break;
                 case PercentileMethod.Nearest:
-                    result = Nearest(percentile);
+                    return SortedArrayStatistics.QuantileCompatible(_data, percentile, QuantileCompatibility.R3);
-                    break;
                 case PercentileMethod.Interpolation:
-                    result = Interpolation(percentile);
+                    return SortedArrayStatistics.QuantileCompatible(_data, percentile, QuantileCompatibility.R5);
-                    break;
                 case PercentileMethod.Excel:
-                    result = Excel(percentile);
+                    return SortedArrayStatistics.QuantileCompatible(_data, percentile, QuantileCompatibility.Excel);
-                    break;
+                default:
+                    return SortedArrayStatistics.Quantile(_data, percentile);
             }
-
-            return result;
         }
 
         /// <summary>
@@ -155,64 +135,5 @@ namespace MathNet.Numerics.Statistics
 
             return percentiles.Select(Compute).ToList();
         }
-
-        /// <summary>
-        /// Computes the percentile using the nearest value.
-        /// </summary>
-        /// <param name="percentile">The percentile.</param>
-        /// <returns>the percentile using the nearest value.</returns>
-        private double Nearest(double percentile)
-        {
-            var n = (int)Math.Round((_data.Count * percentile) + 0.5, 0);
-            return _data[n - 1];
-        }
-
-        /// <summary>
-        /// Computes the percentile using Excel's method.
-        /// </summary>
-        /// <param name="percentile">The percentile.</param>
-        /// <returns>the percentile using Excel's method.</returns>
-        private double Excel(double percentile)
-        {
-            var tmp = 1 + (percentile * (_data.Count - 1.0));
-            var k = (int)tmp;
-            var d = tmp - k;
-
-            return _data[k - 1] + (d * (_data[k] - _data[k - 1]));
-        }
-
-        /// <summary>
-        /// Computes the percentile using interpolation.
-        /// </summary>
-        /// <param name="percentile">The percentile.</param>
-        /// <returns>the percentile using the interpolation.</returns>
-        private double Interpolation(double percentile)
-        {
-            var k = (int)(_data.Count * percentile);
-            var pk = (k - 0.5) / _data.Count;
-            if(k == 0)
-                return _data[0];
-
-            return _data[k - 1] + (_data.Count * (percentile - pk) * (_data[k] - _data[k - 1]));
-        }
-
-        /// <summary>
-        /// Computes the percentile using NIST's method.
-        /// </summary>
-        /// <param name="percentile">The percentile.</param>
-        /// <returns>the percentile using NIST's method.</returns>
-        private double Nist(double percentile)
-        {
-            var tmp = percentile * (_data.Count + 1.0);
-            var k = (int)tmp;
-            if(k == 0)
-                return _data[0];
-            if(k == _data.Count)
-                return _data[k - 1];
-
-            var d = tmp - k;
-
-            return _data[k - 1] + (d * (_data[k] - _data[k - 1]));
-        }
     }
 }

src/Numerics/Statistics/SortedArrayStatistics.cs

@@ -0,0 +1,137 @@
+// <copyright file="SortedArrayStatistics.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>
+
+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
+    }
+
+    public static class SortedArrayStatistics
+    {
+        const double Third = 1d / 3d;
+        const double Half = 1d / 2d;
+
+        /// <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 Quantile(double[] data, double tau)
+        {
+            if (tau < 0d || tau > 1d || data == null || data.Length == 0) return double.NaN;
+            if (tau == 0d || data.Length == 1) return data[0];
+            if (tau == 1d) return data[data.Length - 1];
+
+            double h = (data.Length + Third)*tau + Third;
+            var hf = (int) h;
+            return data[hf - 1] + (h - hf)*(data[hf] - data[hf - 1]);
+        }
+
+        public static double QuantileCompatible(double[] data, double tau, QuantileCompatibility compatibility)
+        {
+            if (tau < 0d || tau > 1d || data == null || 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)
+            {
+                case QuantileCompatibility.R1:
+                case QuantileCompatibility.SAS3:
+                    {
+                        double h = data.Length*tau + Half;
+                        return data[(int) Math.Ceiling(h - Half) - 1];
+                    }
+                case QuantileCompatibility.R2:
+                case QuantileCompatibility.SAS5:
+                    {
+                        double h = data.Length * tau + Half;
+                        return (data[(int) Math.Ceiling(h - Half) - 1] + data[(int) (h + Half) - 1])*Half;
+                    }
+                case QuantileCompatibility.R3:
+                case QuantileCompatibility.SAS2:
+                case QuantileCompatibility.Nearest:
+                    {
+                        double h = data.Length*tau;
+                        return data[(int) Math.Round(h) - 1];
+                    }
+                case QuantileCompatibility.R4:
+                case QuantileCompatibility.SAS1:
+                    {
+                        double h = data.Length*tau;
+                        var hf = (int)h;
+                        return data[hf - 1] + (h - hf) * (data[hf] - data[hf - 1]);
+                    }
+                case QuantileCompatibility.R5:
+                    {
+                        double h = data.Length*tau + Half;
+                        var hf = (int)h;
+                        return data[hf - 1] + (h - hf) * (data[hf] - data[hf - 1]);
+                    }
+                case QuantileCompatibility.R6:
+                case QuantileCompatibility.SAS4:
+                case QuantileCompatibility.Nist:
+                    {
+                        double h = (data.Length + 1)*tau;
+                        var hf = (int)h;
+                        return data[hf - 1] + (h - hf) * (data[hf] - data[hf - 1]);
+                    }
+                case QuantileCompatibility.R7:
+                case QuantileCompatibility.Excel:
+                    {
+                        double h = (data.Length - 1)*tau + 1d;
+                        var hf = (int)h;
+                        return data[hf - 1] + (h - hf) * (data[hf] - data[hf - 1]);
+                    }
+                case QuantileCompatibility.R8:
+                case QuantileCompatibility.Default:
+                    {
+                        double h = (data.Length + Third) * tau + Third;
+                        var hf = (int)h;
+                        return data[hf - 1] + (h - hf) * (data[hf] - data[hf - 1]);
+                    }
+                case QuantileCompatibility.R9:
+                    {
+                        double h = (data.Length + 1d/4d) * tau + 3d/8d;
+                        var hf = (int)h;
+                        return data[hf - 1] + (h - hf) * (data[hf] - data[hf - 1]);
+                    }
+                default:
+                    throw new NotSupportedException();
+            }
+        }
+    }
+}

src/Portable/Portable.csproj

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