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