Initial version of Moserware.Skills TrueSkill calculator to go along with my Computing Your Skill blog post

This commit is contained in:
Jeff Moser
2010-03-18 07:39:48 -04:00
commit cb46631ff8
77 changed files with 6172 additions and 0 deletions

View File

@ -0,0 +1,240 @@
using System;
namespace Moserware.Numerics
{
public class GaussianDistribution
{
// Intentionally, we're not going to derive related things, but set them all at once
// to get around some NaN issues
private GaussianDistribution()
{
}
public GaussianDistribution(double mean, double standardDeviation)
{
Mean = mean;
StandardDeviation = standardDeviation;
Variance = Square(StandardDeviation);
Precision = 1.0/Variance;
PrecisionMean = Precision*Mean;
}
public double Mean { get; private set; }
public double StandardDeviation { get; private set; }
// Precision and PrecisionMean are used because they make multiplying and dividing simpler
// (the the accompanying math paper for more details)
public double Precision { get; private set; }
public double PrecisionMean { get; private set; }
private double Variance { get; set; }
public double NormalizationConstant
{
get
{
// Great derivation of this is at http://www.astro.psu.edu/~mce/A451_2/A451/downloads/notes0.pdf
return 1.0/(Math.Sqrt(2*Math.PI)*StandardDeviation);
}
}
public GaussianDistribution Clone()
{
var result = new GaussianDistribution();
result.Mean = Mean;
result.StandardDeviation = StandardDeviation;
result.Variance = Variance;
result.Precision = Precision;
result.PrecisionMean = PrecisionMean;
return result;
}
public static GaussianDistribution FromPrecisionMean(double precisionMean, double precision)
{
var gaussianDistribution = new GaussianDistribution();
gaussianDistribution.Precision = precision;
gaussianDistribution.PrecisionMean = precisionMean;
gaussianDistribution.Variance = 1.0/precision;
gaussianDistribution.StandardDeviation = Math.Sqrt(gaussianDistribution.Variance);
gaussianDistribution.Mean = gaussianDistribution.PrecisionMean/gaussianDistribution.Precision;
return gaussianDistribution;
}
// Although we could use equations from // For details, see http://www.tina-vision.net/tina-knoppix/tina-memo/2003-003.pdf
// for multiplication, the precision mean ones are easier to write :)
public static GaussianDistribution operator *(GaussianDistribution left, GaussianDistribution right)
{
return FromPrecisionMean(left.PrecisionMean + right.PrecisionMean, left.Precision + right.Precision);
}
/// Computes the absolute difference between two Gaussians
public static double AbsoluteDifference(GaussianDistribution left, GaussianDistribution right)
{
return Math.Max(
Math.Abs(left.PrecisionMean - right.PrecisionMean),
Math.Sqrt(Math.Abs(left.Precision - right.Precision)));
}
/// Computes the absolute difference between two Gaussians
public static double operator -(GaussianDistribution left, GaussianDistribution right)
{
return AbsoluteDifference(left, right);
}
public static double LogProductNormalization(GaussianDistribution left, GaussianDistribution right)
{
if ((left.Precision == 0) || (right.Precision == 0))
{
return 0;
}
double varianceSum = left.Variance + right.Variance;
double meanDifference = left.Mean - right.Mean;
double logSqrt2Pi = Math.Log(Math.Sqrt(2*Math.PI));
return -logSqrt2Pi - (Math.Log(varianceSum)/2.0) - (Square(meanDifference)/(2.0*varianceSum));
}
public static GaussianDistribution operator /(GaussianDistribution numerator, GaussianDistribution denominator)
{
return FromPrecisionMean(numerator.PrecisionMean - denominator.PrecisionMean,
numerator.Precision - denominator.Precision);
}
public static double LogRatioNormalization(GaussianDistribution numerator, GaussianDistribution denominator)
{
if ((numerator.Precision == 0) || (denominator.Precision == 0))
{
return 0;
}
double varianceDifference = denominator.Variance - numerator.Variance;
double meanDifference = numerator.Mean - denominator.Mean;
double logSqrt2Pi = Math.Log(Math.Sqrt(2*Math.PI));
return Math.Log(denominator.Variance) + logSqrt2Pi - Math.Log(varianceDifference)/2.0 +
Square(meanDifference)/(2*varianceDifference);
}
private static double Square(double x)
{
return x*x;
}
public static double At(double x)
{
return At(x, 0, 1);
}
public static double At(double x, double mean, double standardDeviation)
{
// See http://mathworld.wolfram.com/NormalDistribution.html
// 1 -(x-mean)^2 / (2*stdDev^2)
// P(x) = ------------------- * e
// stdDev * sqrt(2*pi)
double multiplier = 1.0/(standardDeviation*Math.Sqrt(2*Math.PI));
double expPart = Math.Exp((-1.0*Math.Pow(x - mean, 2.0))/(2*(standardDeviation*standardDeviation)));
double result = multiplier*expPart;
return result;
}
public static double CumulativeTo(double x, double mean, double standardDeviation)
{
double invsqrt2 = -0.707106781186547524400844362104;
double result = ErrorFunctionCumulativeTo(invsqrt2*x);
return 0.5*result;
}
public static double CumulativeTo(double x)
{
return CumulativeTo(x, 0, 1);
}
private static double ErrorFunctionCumulativeTo(double x)
{
// Derived from page 265 of Numerical Recipes 3rd Edition
double z = Math.Abs(x);
double t = 2.0/(2.0 + z);
double ty = 4*t - 2;
double[] coefficients = {
-1.3026537197817094, 6.4196979235649026e-1,
1.9476473204185836e-2, -9.561514786808631e-3, -9.46595344482036e-4,
3.66839497852761e-4, 4.2523324806907e-5, -2.0278578112534e-5,
-1.624290004647e-6, 1.303655835580e-6, 1.5626441722e-8, -8.5238095915e-8,
6.529054439e-9, 5.059343495e-9, -9.91364156e-10, -2.27365122e-10,
9.6467911e-11, 2.394038e-12, -6.886027e-12, 8.94487e-13, 3.13092e-13,
-1.12708e-13, 3.81e-16, 7.106e-15, -1.523e-15, -9.4e-17, 1.21e-16, -2.8e-17
};
int ncof = coefficients.Length;
double d = 0.0;
double dd = 0.0;
for (int j = ncof - 1; j > 0; j--)
{
double tmp = d;
d = ty*d - dd + coefficients[j];
dd = tmp;
}
double ans = t*Math.Exp(-z*z + 0.5*(coefficients[0] + ty*d) - dd);
return x >= 0.0 ? ans : (2.0 - ans);
}
private static double InverseErrorFunctionCumulativeTo(double p)
{
// From page 265 of numerical recipes
if (p >= 2.0)
{
return -100;
}
if (p <= 0.0)
{
return 100;
}
double pp = (p < 1.0) ? p : 2 - p;
double t = Math.Sqrt(-2*Math.Log(pp/2.0)); // Initial guess
double x = -0.70711*((2.30753 + t*0.27061)/(1.0 + t*(0.99229 + t*0.04481)) - t);
for (int j = 0; j < 2; j++)
{
double err = ErrorFunctionCumulativeTo(x) - pp;
x += err/(1.12837916709551257*Math.Exp(-(x*x)) - x*err); // Halley
}
return p < 1.0 ? x : -x;
}
public static double InverseCumulativeTo(double x, double mean, double standardDeviation)
{
// From numerical recipes, page 320
return mean - Math.Sqrt(2)*standardDeviation*InverseErrorFunctionCumulativeTo(2*x);
}
public static double InverseCumulativeTo(double x)
{
return InverseCumulativeTo(x, 0, 1);
}
public override string ToString()
{
// Debug help
return String.Format("μ={0:0.0000}, σ={1:0.0000}",
Mean,
StandardDeviation);
}
}
}

520
Skills/Numerics/Matrix.cs Normal file
View File

@ -0,0 +1,520 @@
using System;
using System.Collections.Generic;
using System.Linq;
namespace Moserware.Numerics
{
/// <summary>
/// Represents an MxN matrix with double precision values.
/// </summary>
internal class Matrix
{
protected double[][] _MatrixRowValues;
// Note: some properties like Determinant, Inverse, etc are properties instead
// of methods to make the syntax look nicer even though this sort of goes against
// Framework Design Guidelines that properties should be "cheap" since it could take
// a long time to compute these properties if the matrices are "big."
protected Matrix()
{
}
public Matrix(int rows, int columns, params double[] allRowValues)
{
Rows = rows;
Columns = columns;
_MatrixRowValues = new double[rows][];
int currentIndex = 0;
for (int currentRow = 0; currentRow < Rows; currentRow++)
{
_MatrixRowValues[currentRow] = new double[Columns];
for (int currentColumn = 0; currentColumn < Columns; currentColumn++)
{
if ((allRowValues != null) && (currentIndex < allRowValues.Length))
{
_MatrixRowValues[currentRow][currentColumn] = allRowValues[currentIndex++];
}
}
}
}
public Matrix(double[][] rowValues)
{
if (!rowValues.All(row => row.Length == rowValues[0].Length))
{
throw new ArgumentException("All rows must be the same length!");
}
Rows = rowValues.Length;
Columns = rowValues[0].Length;
_MatrixRowValues = rowValues;
}
protected Matrix(int rows, int columns, double[][] matrixRowValues)
{
Rows = rows;
Columns = columns;
_MatrixRowValues = matrixRowValues;
}
public Matrix(int rows, int columns, IEnumerable<IEnumerable<double>> columnValues)
: this(rows, columns)
{
int columnIndex = 0;
foreach (var currentColumn in columnValues)
{
int rowIndex = 0;
foreach (double currentColumnValue in currentColumn)
{
_MatrixRowValues[rowIndex++][columnIndex] = currentColumnValue;
}
columnIndex++;
}
}
public int Rows { get; protected set; }
public int Columns { get; protected set; }
public double this[int row, int column]
{
get { return _MatrixRowValues[row][column]; }
}
public Matrix Transpose
{
get
{
// Just flip everything
var transposeMatrix = new double[Columns][];
for (int currentRowTransposeMatrix = 0;
currentRowTransposeMatrix < Columns;
currentRowTransposeMatrix++)
{
var transposeMatrixCurrentRowColumnValues = new double[Rows];
transposeMatrix[currentRowTransposeMatrix] = transposeMatrixCurrentRowColumnValues;
for (int currentColumnTransposeMatrix = 0;
currentColumnTransposeMatrix < Rows;
currentColumnTransposeMatrix++)
{
transposeMatrixCurrentRowColumnValues[currentColumnTransposeMatrix] =
_MatrixRowValues[currentColumnTransposeMatrix][currentRowTransposeMatrix];
}
}
return new Matrix(Columns, Rows, transposeMatrix);
}
}
private bool IsSquare
{
get { return (Rows == Columns) && Rows > 0; }
}
public double Determinant
{
get
{
// Basic argument checking
if (!IsSquare)
{
throw new NotSupportedException("Matrix must be square!");
}
if (Rows == 1)
{
// Really happy path :)
return _MatrixRowValues[0][0];
}
if (Rows == 2)
{
// Happy path!
// Given:
// | a b |
// | c d |
// The determinant is ad - bc
double a = _MatrixRowValues[0][0];
double b = _MatrixRowValues[0][1];
double c = _MatrixRowValues[1][0];
double d = _MatrixRowValues[1][1];
return a*d - b*c;
}
// I use the Laplace expansion here since it's straightforward to implement.
// It's O(n^2) and my implementation is especially poor performing, but the
// core idea is there. Perhaps I should replace it with a better algorithm
// later.
// See http://en.wikipedia.org/wiki/Laplace_expansion for details
double result = 0.0;
// I expand along the first row
for (int currentColumn = 0; currentColumn < Columns; currentColumn++)
{
double firstRowColValue = _MatrixRowValues[0][currentColumn];
double cofactor = GetCofactor(0, currentColumn);
double itemToAdd = firstRowColValue*cofactor;
result += itemToAdd;
}
return result;
}
}
public Matrix Adjugate
{
get
{
if (!IsSquare)
{
throw new ArgumentException("Matrix must be square!");
}
// See http://en.wikipedia.org/wiki/Adjugate_matrix
if (Rows == 2)
{
// Happy path!
// Adjugate of:
// | a b |
// | c d |
// is
// | d -b |
// | -c a |
double a = _MatrixRowValues[0][0];
double b = _MatrixRowValues[0][1];
double c = _MatrixRowValues[1][0];
double d = _MatrixRowValues[1][1];
return new SquareMatrix(d, -b,
-c, a);
}
// The idea is that it's the transpose of the cofactors
var result = new double[Columns][];
for (int currentColumn = 0; currentColumn < Columns; currentColumn++)
{
result[currentColumn] = new double[Rows];
for (int currentRow = 0; currentRow < Rows; currentRow++)
{
result[currentColumn][currentRow] = GetCofactor(currentRow, currentColumn);
}
}
return new Matrix(result);
}
}
public Matrix Inverse
{
get
{
if ((Rows == 1) && (Columns == 1))
{
return new SquareMatrix(1.0/_MatrixRowValues[0][0]);
}
// Take the simple approach:
// http://en.wikipedia.org/wiki/Cramer%27s_rule#Finding_inverse_matrix
return (1.0/Determinant)*Adjugate;
}
}
public static Matrix operator *(double scalarValue, Matrix matrix)
{
int rows = matrix.Rows;
int columns = matrix.Columns;
var newValues = new double[rows][];
for (int currentRow = 0; currentRow < rows; currentRow++)
{
var newRowColumnValues = new double[columns];
newValues[currentRow] = newRowColumnValues;
for (int currentColumn = 0; currentColumn < columns; currentColumn++)
{
newRowColumnValues[currentColumn] = scalarValue*matrix._MatrixRowValues[currentRow][currentColumn];
}
}
return new Matrix(rows, columns, newValues);
}
public static Matrix operator +(Matrix left, Matrix right)
{
if ((left.Rows != right.Rows) || (left.Columns != right.Columns))
{
throw new ArgumentException("Matrices must be of the same size");
}
// simple addition of each item
var resultMatrix = new double[left.Rows][];
for (int currentRow = 0; currentRow < left.Rows; currentRow++)
{
var rowColumnValues = new double[right.Columns];
resultMatrix[currentRow] = rowColumnValues;
for (int currentColumn = 0; currentColumn < right.Columns; currentColumn++)
{
rowColumnValues[currentColumn] = left._MatrixRowValues[currentRow][currentColumn]
+
right._MatrixRowValues[currentRow][currentColumn];
}
}
return new Matrix(left.Rows, right.Columns, resultMatrix);
}
public static Matrix operator *(Matrix left, Matrix right)
{
// Just your standard matrix multiplication.
// See http://en.wikipedia.org/wiki/Matrix_multiplication for details
if (left.Columns != right.Rows)
{
throw new ArgumentException("The width of the left matrix must match the height of the right matrix",
"right");
}
int resultRows = left.Rows;
int resultColumns = right.Columns;
var resultMatrix = new double[resultRows][];
for (int currentRow = 0; currentRow < resultRows; currentRow++)
{
resultMatrix[currentRow] = new double[resultColumns];
for (int currentColumn = 0; currentColumn < resultColumns; currentColumn++)
{
double productValue = 0;
for (int vectorIndex = 0; vectorIndex < left.Columns; vectorIndex++)
{
double leftValue = left._MatrixRowValues[currentRow][vectorIndex];
double rightValue = right._MatrixRowValues[vectorIndex][currentColumn];
double vectorIndexProduct = leftValue*rightValue;
productValue += vectorIndexProduct;
}
resultMatrix[currentRow][currentColumn] = productValue;
}
}
return new Matrix(resultRows, resultColumns, resultMatrix);
}
private Matrix GetMinorMatrix(int rowToRemove, int columnToRemove)
{
// See http://en.wikipedia.org/wiki/Minor_(linear_algebra)
// I'm going to use a horribly naïve algorithm... because I can :)
var result = new double[Rows - 1][];
int resultRow = 0;
for (int currentRow = 0; currentRow < Rows; currentRow++)
{
if (currentRow == rowToRemove)
{
continue;
}
result[resultRow] = new double[Columns - 1];
int resultColumn = 0;
for (int currentColumn = 0; currentColumn < Columns; currentColumn++)
{
if (currentColumn == columnToRemove)
{
continue;
}
result[resultRow][resultColumn] = _MatrixRowValues[currentRow][currentColumn];
resultColumn++;
}
resultRow++;
}
return new Matrix(Rows - 1, Columns - 1, result);
}
private double GetCofactor(int rowToRemove, int columnToRemove)
{
// See http://en.wikipedia.org/wiki/Cofactor_(linear_algebra) for details
// REVIEW: should things be reversed since I'm 0 indexed?
int sum = rowToRemove + columnToRemove;
bool isEven = (sum%2 == 0);
if (isEven)
{
return GetMinorMatrix(rowToRemove, columnToRemove).Determinant;
}
else
{
return -1.0*GetMinorMatrix(rowToRemove, columnToRemove).Determinant;
}
}
// Equality stuff
// See http://msdn.microsoft.com/en-us/library/ms173147.aspx
public static bool operator ==(Matrix a, Matrix b)
{
// If both are null, or both are same instance, return true.
if (ReferenceEquals(a, b))
{
return true;
}
// If one is null, but not both, return false.
if (((object) a == null) || ((object) b == null))
{
return false;
}
if ((a.Rows != b.Rows) || (a.Columns != b.Columns))
{
return false;
}
const double errorTolerance = 0.0000000000001;
for (int currentRow = 0; currentRow < a.Rows; currentRow++)
{
for (int currentColumn = 0; currentColumn < a.Columns; currentColumn++)
{
double delta =
Math.Abs(a._MatrixRowValues[currentRow][currentColumn] -
b._MatrixRowValues[currentRow][currentColumn]);
if (delta > errorTolerance)
{
return false;
}
}
}
return true;
}
public static bool operator !=(Matrix a, Matrix b)
{
return !(a == b);
}
public override int GetHashCode()
{
double result = Rows;
result += 2*Columns;
unchecked
{
for (int currentRow = 0; currentRow < Rows; currentRow++)
{
bool eventRow = (currentRow%2) == 0;
double multiplier = eventRow ? 1.0 : 2.0;
for (int currentColumn = 0; currentColumn < Columns; currentColumn++)
{
result += multiplier*_MatrixRowValues[currentRow][currentColumn];
}
}
}
// Ok, now convert that double to an int
byte[] resultBytes = BitConverter.GetBytes(result);
var finalBytes = new byte[4];
for (int i = 0; i < 4; i++)
{
finalBytes[i] = (byte) (resultBytes[i] ^ resultBytes[i + 4]);
}
int hashCode = BitConverter.ToInt32(finalBytes, 0);
return hashCode;
}
public override bool Equals(object obj)
{
var other = obj as Matrix;
if (other == null)
{
return base.Equals(obj);
}
return this == other;
}
}
internal class DiagonalMatrix : Matrix
{
public DiagonalMatrix(IList<double> diagonalValues)
: base(diagonalValues.Count, diagonalValues.Count)
{
for (int i = 0; i < diagonalValues.Count; i++)
{
_MatrixRowValues[i][i] = diagonalValues[i];
}
}
}
internal class Vector : Matrix
{
public Vector(IList<double> vectorValues)
: base(vectorValues.Count, 1, new IEnumerable<double>[] {vectorValues})
{
}
}
internal class SquareMatrix : Matrix
{
public SquareMatrix(params double[] allValues)
{
Rows = (int) Math.Sqrt(allValues.Length);
Columns = Rows;
int allValuesIndex = 0;
_MatrixRowValues = new double[Rows][];
for (int currentRow = 0; currentRow < Rows; currentRow++)
{
var currentRowValues = new double[Columns];
_MatrixRowValues[currentRow] = currentRowValues;
for (int currentColumn = 0; currentColumn < Columns; currentColumn++)
{
currentRowValues[currentColumn] = allValues[allValuesIndex++];
}
}
}
}
internal class IdentityMatrix : DiagonalMatrix
{
public IdentityMatrix(int rows)
: base(CreateDiagonal(rows))
{
}
private static double[] CreateDiagonal(int rows)
{
var result = new double[rows];
for (int i = 0; i < rows; i++)
{
result[i] = 1.0;
}
return result;
}
}
}

49
Skills/Numerics/Range.cs Normal file
View File

@ -0,0 +1,49 @@
using System;
namespace Moserware.Skills.Numerics
{
// The whole purpose of this class is to make the code for the SkillCalculator(s)
// look a little cleaner
public abstract class Range<T> where T : Range<T>, new()
{
private static readonly T _Instance = new T();
protected Range(int min, int max)
{
if (min > max)
{
throw new ArgumentOutOfRangeException();
}
Min = min;
Max = max;
}
public int Min { get; private set; }
public int Max { get; private set; }
protected abstract T Create(int min, int max);
// REVIEW: It's probably bad form to have access statics via a derived class, but the syntax looks better :-)
public static T Inclusive(int min, int max)
{
return _Instance.Create(min, max);
}
public static T Exactly(int value)
{
return _Instance.Create(value, value);
}
public static T AtLeast(int minimumValue)
{
return _Instance.Create(minimumValue, int.MaxValue);
}
public bool IsInRange(int value)
{
return (Min <= value) && (value <= Max);
}
}
}