tsai
/
mathnet-numerics
1
0
Fork 3
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

Statistics: Ranks

pull/194/head
Christoph Ruegg 13 years ago
parent
3aca1f4bf7
commit
0b5f954739
9 changed files with 638 additions and 274 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 1
      src/Numerics/Numerics.csproj
  2. 265
      src/Numerics/Sorting.cs
  3. 210
      src/Numerics/Statistics/ArrayStatistics.cs
  4. 37
      src/Numerics/Statistics/Correlation.cs
  5. 49
      src/Numerics/Statistics/RankDefinition.cs
  6. 196
      src/Numerics/Statistics/SortedArrayStatistics.cs
  7. 26
      src/Numerics/Statistics/Statistics.cs
  8. 7
      src/UnitTests/ExcelTests.cs
  9. 121
      src/UnitTests/StatisticsTests/StatisticsTests.cs

1
src/Numerics/Numerics.csproj
View File

@ -202,6 +202,7 @@
<Compile Include="SpecialFunctions\Logistic.cs" />
<Compile Include="Statistics\ArrayStatistics.cs" />
<Compile Include="Statistics\QuantileDefinition.cs" />
<Compile Include="Statistics\RankDefinition.cs" />
<Compile Include="Statistics\StreamingStatistics.cs" />
<Compile Include="Statistics\SortedArrayStatistics.cs" />
<Compile Include="LinearAlgebra\Storage\SparseVectorStorage.cs" />

265
src/Numerics/Sorting.cs
View File

@ -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.
// </copyright>
using System;
using System.Collections.Generic;
namespace MathNet.Numerics
{
using System;
using System.Collections.Generic;
/// <summary>
/// Sorting algorithms for single, tuple and triple lists.
/// </summary>
public static class Sorting
{
/// <summary>
/// Sort a list of keys, in place using the quick sort algorithm.
/// </summary>
/// <typeparam name="T">The type of elements stored in the list.</typeparam>
/// <param name="keys">List to sort.</param>
public static void Sort<T>(IList<T> keys)
{
Sort(keys, Comparer<T>.Default);
}
/// <summary>
/// Sort a list of keys and items with respect to the keys, in place using the quick sort algorithm.
/// </summary>
/// <typeparam name="TKey">The type of elements stored in the key list.</typeparam>
/// <typeparam name="TItem">The type of elements stored in the item list.</typeparam>
/// <param name="keys">List to sort.</param>
/// <param name="items">List to permute the same way as the key list.</param>
public static void Sort<TKey, TItem>(IList<TKey> keys, IList<TItem> items)
{
Sort(keys, items, Comparer<TKey>.Default);
}
/// <summary>
/// Sort a list of keys, items1 and items2 with respect to the keys, in place using the quick sort algorithm.
/// </summary>
/// <typeparam name="TKey">The type of elements stored in the key list.</typeparam>
/// <typeparam name="TItem1">The type of elements stored in the first item list.</typeparam>
/// <typeparam name="TItem2">The type of elements stored in the second item list.</typeparam>
/// <param name="keys">List to sort.</param>
/// <param name="items1">First list to permute the same way as the key list.</param>
/// <param name="items2">Second list to permute the same way as the key list.</param>
public static void Sort<TKey, TItem1, TItem2>(IList<TKey> keys, IList<TItem1> items1, IList<TItem2> items2)
{
Sort(keys, items1, items2, Comparer<TKey>.Default);
}
/// <summary>
/// Sort a range of a list of keys, in place using the quick sort algorithm.
/// </summary>
@ -92,16 +56,11 @@ namespace MathNet.Numerics
/// <typeparam name="T">The type of elements in the key list.</typeparam>
/// <param name="keys">List to sort.</param>
/// <param name="comparer">Comparison, defining the sort order.</param>
public static void Sort<T>(IList<T> keys, IComparer<T> comparer)
public static void Sort<T>(IList<T> keys, IComparer<T> comparer = null)
{
if (null == keys)
{
throw new ArgumentNullException("keys");
}
if (null == comparer)
{
throw new ArgumentNullException("comparer");
comparer = Comparer<T>.Default;
}
// basic cases
@ -148,21 +107,11 @@ namespace MathNet.Numerics
/// <param name="keys">List to sort.</param>
/// <param name="items">List to permute the same way as the key list.</param>
/// <param name="comparer">Comparison, defining the sort order.</param>
public static void Sort<TKey, TItem>(IList<TKey> keys, IList<TItem> items, IComparer<TKey> comparer)
public static void Sort<TKey, TItem>(IList<TKey> keys, IList<TItem> items, IComparer<TKey> 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<TKey>.Default;
}
#if !PORTABLE
@ -190,31 +139,40 @@ namespace MathNet.Numerics
/// <param name="items1">First list to permute the same way as the key list.</param>
/// <param name="items2">Second list to permute the same way as the key list.</param>
/// <param name="comparer">Comparison, defining the sort order.</param>
public static void Sort<TKey, TItem1, TItem2>(
IList<TKey> keys, IList<TItem1> items1, IList<TItem2> items2, IComparer<TKey> comparer)
public static void Sort<TKey, TItem1, TItem2>(IList<TKey> keys, IList<TItem1> items1, IList<TItem2> items2, IComparer<TKey> comparer = null)
{
if (null == keys)
if (null == comparer)
{
throw new ArgumentNullException("keys");
comparer = Comparer<TKey>.Default;
}
if (null == items1)
{
throw new ArgumentNullException("items1");
}
// local sort implementation
QuickSort(keys, items1, items2, comparer, 0, keys.Count - 1);
}
if (null == items2)
/// <summary>
/// Sort a list of keys and items with respect to the keys, in place using the quick sort algorithm.
/// </summary>
/// <typeparam name="T1">The type of elements in the primary list.</typeparam>
/// <typeparam name="T2">The type of elements in the secondary list.</typeparam>
/// <param name="primary">List to sort.</param>
/// <param name="secondary">List to to sort on duplicate primaty items, and permute the same way as the key list.</param>
/// <param name="primaryComparer">Comparison, defining the primary sort order.</param>
/// <param name="secondaryComparer">Comparison, defining the secondary sort order.</param>
public static void SortAll<T1, T2>(IList<T1> primary, IList<T2> secondary, IComparer<T1> primaryComparer = null, IComparer<T2> secondaryComparer = null)
{
if (null == primaryComparer)
{
throw new ArgumentNullException("items2");
primaryComparer = Comparer<T1>.Default;
}
if (null == comparer)
if (null == secondaryComparer)
{
throw new ArgumentNullException("comparer");
secondaryComparer = Comparer<T2>.Default;
}
// local sort implementation
QuickSort(keys, items1, items2, comparer, 0, keys.Count - 1);
QuickSortAll(primary, secondary, primaryComparer, secondaryComparer, 0, primary.Count - 1);
}
/// <summary>
@ -225,18 +183,8 @@ namespace MathNet.Numerics
/// <param name="index">The zero-based starting index of the range to sort.</param>
/// <param name="count">The length of the range to sort.</param>
/// <param name="comparer">Comparison, defining the sort order.</param>
public static void Sort<T>(IList<T> keys, int index, int count, IComparer<T> comparer)
public static void Sort<T>(IList<T> keys, int index, int count, IComparer<T> 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<T>.Default;
}
// basic cases
if (count <= 1)
{
@ -291,11 +244,7 @@ namespace MathNet.Numerics
/// <param name="comparer">The method with which to compare two elements of the quick sort.</param>
/// <param name="left">The left boundary of the quick sort.</param>
/// <param name="right">The right boundary of the quick sort.</param>
private static void QuickSort<T>(
IList<T> keys,
IComparer<T> comparer,
int left,
int right)
static void QuickSort<T>(IList<T> keys, IComparer<T> 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);
}
/// <summary>
@ -383,12 +330,7 @@ namespace MathNet.Numerics
/// <param name="comparer">The method with which to compare two elements of the quick sort.</param>
/// <param name="left">The left boundary of the quick sort.</param>
/// <param name="right">The right boundary of the quick sort.</param>
private static void QuickSort<T, TItems>(
IList<T> keys,
IList<TItems> items,
IComparer<T> comparer,
int left,
int right)
static void QuickSort<T, TItems>(IList<T> keys, IList<TItems> items, IComparer<T> 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);
}
/// <summary>
@ -482,13 +422,10 @@ namespace MathNet.Numerics
/// <param name="comparer">The method with which to compare two elements of the quick sort.</param>
/// <param name="left">The left boundary of the quick sort.</param>
/// <param name="right">The right boundary of the quick sort.</param>
private static void QuickSort<T, TItems1, TItems2>(
IList<T> keys,
IList<TItems1> items1,
IList<TItems2> items2,
static void QuickSort<T, TItems1, TItems2>(
IList<T> keys, IList<TItems1> items1, IList<TItems2> items2,
IComparer<T> 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);
}
/// <summary>
/// Recursive implementation for an in place quick sort on the primary and then by the secondary list while reordering one secondary list accordingly.
/// </summary>
/// <typeparam name="T1">The type of the primary list.</typeparam>
/// <typeparam name="T2">The type of the secondary list.</typeparam>
/// <param name="primary">The list which is sorted using quick sort.</param>
/// <param name="secondary">The list which is sorted secondarily (on primary duplicates) and automatically reordered accordingly.</param>
/// <param name="primaryComparer">The method with which to compare two elements of the primary list.</param>
/// <param name="secondaryComparer">The method with which to compare two elements of the secondary list.</param>
/// <param name="left">The left boundary of the quick sort.</param>
/// <param name="right">The right boundary of the quick sort.</param>
static void QuickSortAll<T1, T2>(
IList<T1> primary, IList<T2> secondary,
IComparer<T1> primaryComparer, IComparer<T2> 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);
}
/// <summary>

210
src/Numerics/Statistics/ArrayStatistics.cs
View File

@ -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
/// <summary>
/// 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
/// <param name="p">Percentile selector, between 0 and 100 (inclusive).</param>
public static double PercentileInplace(double[] data, int p)
{
return QuantileInplace(data, p / 100d);
return QuantileInplace(data, p/100d);
}
/// <summary>
@ -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);
}
/// <summary>
@ -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;
}
}
/// <summary>
/// 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.
/// </summary>
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;
}
}
}
}
}

37
src/Numerics/Statistics/Correlation.cs
View File

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

49
src/Numerics/Statistics/RankDefinition.cs
View File

@ -0,0 +1,49 @@
// <copyright file="RankDefinition.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-2013 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use,
// copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following
// conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>
namespace MathNet.Numerics.Statistics
{
public enum RankDefinition
{
/// <summary>Replace ties with their mean (non-integer ranks). Default.</summary>
Average = 1,
Default = 1,
/// <summary>Replace ties with their minimum (typical sports ranking).</summary>
Min = 2,
Sports = 2,
/// <summary>Replace ties with their maximum.</summary>
Max = 3,
/// <summary>Permutation with increasing values at each index of ties.</summary>
First = 4
}
}

196
src/Numerics/Statistics/SortedArrayStatistics.cs
View File

@ -133,8 +133,8 @@ namespace MathNet.Numerics.Statistics
/// <param name="data">Sample array, must be sorted ascendingly.</param>
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] };
}
/// <summary>
@ -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]);
}
/// <summary>
@ -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();
}
}
/// <summary>
/// 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.
/// </summary>
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;
}
}
}
}

26
src/Numerics/Statistics/Statistics.cs
View File

@ -634,5 +634,31 @@ namespace MathNet.Numerics.Statistics
Array.Sort(array);
return order => SortedArrayStatistics.OrderStatistic(array, order);
}
/// <summary>
/// Evaluates the rank of each entry of the provided samples.
/// The rank definition can be specificed to be compatible
/// with an existing system.
/// </summary>
/// <param name="data">The data sample sequence.</param>
/// <param name="definition">Rank definition, to choose how ties should be handled.</param>
public static double[] Ranks(this IEnumerable<double> data, RankDefinition definition = RankDefinition.Default)
{
var array = data.ToArray();
return ArrayStatistics.RanksInplace(array, definition);
}
/// <summary>
/// Evaluates the rank of each entry of the provided samples.
/// The rank definition can be specificed to be compatible
/// with an existing system.
/// </summary>
/// <param name="data">The data sample sequence.</param>
/// <param name="definition">Rank definition, to choose how ties should be handled.</param>
public static double[] Ranks(this IEnumerable<double?> data, RankDefinition definition = RankDefinition.Default)
{
var array = data.Where(d => d.HasValue).Select(d => d.Value).ToArray();
return ArrayStatistics.RanksInplace(array, definition);
}
}
}

7
src/UnitTests/ExcelTests.cs
View File

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

121
src/UnitTests/StatisticsTests/StatisticsTests.cs
View File

@ -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

Write Preview
Loading…
Cancel
Save