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