forked from tsai/mathnet-numerics
Christoph Ruegg
13 years ago
8 changed files with 430 additions and 64 deletions
@ -0,0 +1,169 @@ |
|||
// <copyright file="MultipleRegression.cs" company="Math.NET">
|
|||
// Math.NET Numerics, part of the Math.NET Project
|
|||
// http://numerics.mathdotnet.com
|
|||
// http://github.com/mathnet/mathnet-numerics
|
|||
// http://mathnetnumerics.codeplex.com
|
|||
//
|
|||
// Copyright (c) 2009-2013 Math.NET
|
|||
//
|
|||
// Permission is hereby granted, free of charge, to any person
|
|||
// obtaining a copy of this software and associated documentation
|
|||
// files (the "Software"), to deal in the Software without
|
|||
// restriction, including without limitation the rights to use,
|
|||
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|||
// copies of the Software, and to permit persons to whom the
|
|||
// Software is furnished to do so, subject to the following
|
|||
// conditions:
|
|||
//
|
|||
// The above copyright notice and this permission notice shall be
|
|||
// included in all copies or substantial portions of the Software.
|
|||
//
|
|||
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|||
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
|||
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
|||
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
|||
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
|||
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|||
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
|||
// OTHER DEALINGS IN THE SOFTWARE.
|
|||
// </copyright>
|
|||
|
|||
using System; |
|||
using System.Collections.Generic; |
|||
using MathNet.Numerics.LinearAlgebra; |
|||
|
|||
namespace MathNet.Numerics.LinearRegression |
|||
{ |
|||
public static class MultipleRegression |
|||
{ |
|||
/// <summary>
|
|||
/// Find the model parameters β such that X*β with predictor X becomes as close to response Y as possible, with least squares residuals.
|
|||
/// Uses the cholesky decomposition of the normal equations.
|
|||
/// </summary>
|
|||
/// <param name="x">Predictor matrix X</param>
|
|||
/// <param name="y">Response vector Y</param>
|
|||
/// <returns>Best fitting vector for model parameters β</returns>
|
|||
public static Vector<T> NormalEquations<T>(Matrix<T> x, Vector<T> y) where T : struct, IEquatable<T>, IFormattable |
|||
{ |
|||
return x.TransposeThisAndMultiply(x).Cholesky().Solve(x.Transpose()*y); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Find the model parameters β such that their linear combination with all predictor-arrays in X become as close to their response in Y as possible, with least squares residuals.
|
|||
/// Uses the cholesky decomposition of the normal equations.
|
|||
/// </summary>
|
|||
/// <param name="x">List of predictor-arrays.</param>
|
|||
/// <param name="y">List of responses</param>
|
|||
/// <param name="intercept">True if an intercept should be added as first artificial perdictor value. Default = false.</param>
|
|||
/// <returns>Best fitting list of model parameters β for each element in the predictor-arrays.</returns>
|
|||
public static T[] NormalEquations<T>(T[][] x, T[] y, bool intercept = false) where T : struct, IEquatable<T>, IFormattable |
|||
{ |
|||
var predictor = Matrix<T>.Build.DenseMatrixOfRowArrays(x); |
|||
if (intercept) |
|||
{ |
|||
predictor = predictor.InsertColumn(0, Vector<T>.Build.DenseVector(predictor.RowCount, Vector<T>.One)); |
|||
} |
|||
var response = Matrix<T>.Build.DenseVector(y); |
|||
return predictor.TransposeThisAndMultiply(predictor).Cholesky().Solve(predictor.Transpose()*response).ToArray(); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Find the model parameters β such that their linear combination with all predictor-arrays in X become as close to their response in Y as possible, with least squares residuals.
|
|||
/// Uses the cholesky decomposition of the normal equations.
|
|||
/// </summary>
|
|||
/// <param name="samples">Sequence of predictor-arrays and their response.</param>
|
|||
/// <param name="intercept">True if an intercept should be added as first artificial perdictor value. Default = false.</param>
|
|||
/// <returns>Best fitting list of model parameters β for each element in the predictor-arrays.</returns>
|
|||
public static T[] NormalEquations<T>(IEnumerable<Tuple<T[], T>> samples, bool intercept = false) where T : struct, IEquatable<T>, IFormattable |
|||
{ |
|||
var xy = samples.UnpackSinglePass(); |
|||
return NormalEquations(xy.Item1, xy.Item2, intercept); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Find the model parameters β such that X*β with predictor X becomes as close to response Y as possible, with least squares residuals.
|
|||
/// Uses an orthogonal decomposition and is therefore more numerically stable than the normal equations but also slower.
|
|||
/// </summary>
|
|||
/// <param name="x">Predictor matrix X</param>
|
|||
/// <param name="y">Response vector Y</param>
|
|||
/// <returns>Best fitting vector for model parameters β</returns>
|
|||
public static Vector<T> QR<T>(Matrix<T> x, Vector<T> y) where T : struct, IEquatable<T>, IFormattable |
|||
{ |
|||
return x.QR().Solve(y); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Find the model parameters β such that their linear combination with all predictor-arrays in X become as close to their response in Y as possible, with least squares residuals.
|
|||
/// Uses an orthogonal decomposition and is therefore more numerically stable than the normal equations but also slower.
|
|||
/// </summary>
|
|||
/// <param name="x">List of predictor-arrays.</param>
|
|||
/// <param name="y">List of responses</param>
|
|||
/// <param name="intercept">True if an intercept should be added as first artificial perdictor value. Default = false.</param>
|
|||
/// <returns>Best fitting list of model parameters β for each element in the predictor-arrays.</returns>
|
|||
public static T[] QR<T>(T[][] x, T[] y, bool intercept = false) where T : struct, IEquatable<T>, IFormattable |
|||
{ |
|||
var predictor = Matrix<T>.Build.DenseMatrixOfRowArrays(x); |
|||
if (intercept) |
|||
{ |
|||
predictor = predictor.InsertColumn(0, Vector<T>.Build.DenseVector(predictor.RowCount, Vector<T>.One)); |
|||
} |
|||
return predictor.QR().Solve(Matrix<T>.Build.DenseVector(y)).ToArray(); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Find the model parameters β such that their linear combination with all predictor-arrays in X become as close to their response in Y as possible, with least squares residuals.
|
|||
/// Uses an orthogonal decomposition and is therefore more numerically stable than the normal equations but also slower.
|
|||
/// </summary>
|
|||
/// <param name="samples">Sequence of predictor-arrays and their response.</param>
|
|||
/// <param name="intercept">True if an intercept should be added as first artificial perdictor value. Default = false.</param>
|
|||
/// <returns>Best fitting list of model parameters β for each element in the predictor-arrays.</returns>
|
|||
public static T[] QR<T>(IEnumerable<Tuple<T[], T>> samples, bool intercept = false) where T : struct, IEquatable<T>, IFormattable |
|||
{ |
|||
var xy = samples.UnpackSinglePass(); |
|||
return QR(xy.Item1, xy.Item2, intercept); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Find the model parameters β such that X*β with predictor X becomes as close to response Y as possible, with least squares residuals.
|
|||
/// Uses a singular value decomposition and is therefore more numerically stable (especially if ill-conditioned) than the normal equations or QR but also slower.
|
|||
/// </summary>
|
|||
/// <param name="x">Predictor matrix X</param>
|
|||
/// <param name="y">Response vector Y</param>
|
|||
/// <returns>Best fitting vector for model parameters β</returns>
|
|||
public static Vector<T> Svd<T>(Matrix<T> x, Vector<T> y) where T : struct, IEquatable<T>, IFormattable |
|||
{ |
|||
return x.Svd().Solve(y); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Find the model parameters β such that their linear combination with all predictor-arrays in X become as close to their response in Y as possible, with least squares residuals.
|
|||
/// Uses a singular value decomposition and is therefore more numerically stable (especially if ill-conditioned) than the normal equations or QR but also slower.
|
|||
/// </summary>
|
|||
/// <param name="x">List of predictor-arrays.</param>
|
|||
/// <param name="y">List of responses</param>
|
|||
/// <param name="intercept">True if an intercept should be added as first artificial perdictor value. Default = false.</param>
|
|||
/// <returns>Best fitting list of model parameters β for each element in the predictor-arrays.</returns>
|
|||
public static T[] Svd<T>(T[][] x, T[] y, bool intercept = false) where T : struct, IEquatable<T>, IFormattable |
|||
{ |
|||
var predictor = Matrix<T>.Build.DenseMatrixOfRowArrays(x); |
|||
if (intercept) |
|||
{ |
|||
predictor = predictor.InsertColumn(0, Vector<T>.Build.DenseVector(predictor.RowCount, Vector<T>.One)); |
|||
} |
|||
return predictor.Svd().Solve(Matrix<T>.Build.DenseVector(y)).ToArray(); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Find the model parameters β such that their linear combination with all predictor-arrays in X become as close to their response in Y as possible, with least squares residuals.
|
|||
/// Uses a singular value decomposition and is therefore more numerically stable (especially if ill-conditioned) than the normal equations or QR but also slower.
|
|||
/// </summary>
|
|||
/// <param name="samples">Sequence of predictor-arrays and their response.</param>
|
|||
/// <param name="intercept">True if an intercept should be added as first artificial perdictor value. Default = false.</param>
|
|||
/// <returns>Best fitting list of model parameters β for each element in the predictor-arrays.</returns>
|
|||
public static T[] Svd<T>(IEnumerable<Tuple<T[], T>> samples, bool intercept = false) where T : struct, IEquatable<T>, IFormattable |
|||
{ |
|||
var xy = samples.UnpackSinglePass(); |
|||
return Svd(xy.Item1, xy.Item2, intercept); |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,88 @@ |
|||
// <copyright file="SimpleRegression.cs" company="Math.NET">
|
|||
// Math.NET Numerics, part of the Math.NET Project
|
|||
// http://numerics.mathdotnet.com
|
|||
// http://github.com/mathnet/mathnet-numerics
|
|||
// http://mathnetnumerics.codeplex.com
|
|||
//
|
|||
// Copyright (c) 2009-2013 Math.NET
|
|||
//
|
|||
// Permission is hereby granted, free of charge, to any person
|
|||
// obtaining a copy of this software and associated documentation
|
|||
// files (the "Software"), to deal in the Software without
|
|||
// restriction, including without limitation the rights to use,
|
|||
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|||
// copies of the Software, and to permit persons to whom the
|
|||
// Software is furnished to do so, subject to the following
|
|||
// conditions:
|
|||
//
|
|||
// The above copyright notice and this permission notice shall be
|
|||
// included in all copies or substantial portions of the Software.
|
|||
//
|
|||
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|||
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
|||
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
|||
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
|||
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
|||
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|||
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
|||
// OTHER DEALINGS IN THE SOFTWARE.
|
|||
// </copyright>
|
|||
|
|||
using System; |
|||
using System.Collections.Generic; |
|||
using MathNet.Numerics.Properties; |
|||
|
|||
namespace MathNet.Numerics.LinearRegression |
|||
{ |
|||
public static class SimpleRegression |
|||
{ |
|||
/// <summary>
|
|||
/// Least-Squares fitting the points (x,y) to a line y : x -> a+b*x,
|
|||
/// returning its best fitting parameters as (a, b) tuple,
|
|||
/// where a is the intercept and b the slope.
|
|||
/// </summary>
|
|||
/// <param name="x">Predictor (independent)</param>
|
|||
/// <param name="y">Response (dependent)</param>
|
|||
public static Tuple<double, double> Fit(double[] x, double[] y) |
|||
{ |
|||
if (x.Length != y.Length) throw new ArgumentException(Resources.ArgumentVectorsSameLength); |
|||
if (x.Length <= 1) throw new ArgumentException(string.Format(Resources.ArrayTooSmall, 2)); |
|||
|
|||
// First Pass: Mean (Less robust but faster than ArrayStatistics.Mean)
|
|||
double mx = 0.0; |
|||
double my = 0.0; |
|||
for (int i = 0; i < x.Length; i++) |
|||
{ |
|||
mx += x[i]; |
|||
my += y[i]; |
|||
} |
|||
mx /= x.Length; |
|||
my /= y.Length; |
|||
|
|||
// Second Pass: Covariance/Variance
|
|||
double covariance = 0.0; |
|||
double variance = 0.0; |
|||
for (int i = 0; i < x.Length; i++) |
|||
{ |
|||
double diff = x[i] - mx; |
|||
covariance += diff*(y[i] - my); |
|||
variance += diff*diff; |
|||
} |
|||
|
|||
var b = covariance/variance; |
|||
return new Tuple<double, double>(my - b*mx, b); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Least-Squares fitting the points (x,y) to a line y : x -> a+b*x,
|
|||
/// returning its best fitting parameters as (a, b) tuple,
|
|||
/// where a is the intercept and b the slope.
|
|||
/// </summary>
|
|||
/// <param name="samples">Predictor-Response samples as tuples</param>
|
|||
public static Tuple<double, double> Fit(IEnumerable<Tuple<double, double>> samples) |
|||
{ |
|||
var xy = samples.UnpackSinglePass(); |
|||
return Fit(xy.Item1, xy.Item2); |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,52 @@ |
|||
// <copyright file="Util.cs" company="Math.NET">
|
|||
// Math.NET Numerics, part of the Math.NET Project
|
|||
// http://numerics.mathdotnet.com
|
|||
// http://github.com/mathnet/mathnet-numerics
|
|||
// http://mathnetnumerics.codeplex.com
|
|||
//
|
|||
// Copyright (c) 2009-2013 Math.NET
|
|||
//
|
|||
// Permission is hereby granted, free of charge, to any person
|
|||
// obtaining a copy of this software and associated documentation
|
|||
// files (the "Software"), to deal in the Software without
|
|||
// restriction, including without limitation the rights to use,
|
|||
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|||
// copies of the Software, and to permit persons to whom the
|
|||
// Software is furnished to do so, subject to the following
|
|||
// conditions:
|
|||
//
|
|||
// The above copyright notice and this permission notice shall be
|
|||
// included in all copies or substantial portions of the Software.
|
|||
//
|
|||
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|||
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
|||
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
|||
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
|||
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
|||
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|||
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
|||
// OTHER DEALINGS IN THE SOFTWARE.
|
|||
// </copyright>
|
|||
|
|||
using System; |
|||
using System.Collections.Generic; |
|||
|
|||
namespace MathNet.Numerics.LinearRegression |
|||
{ |
|||
internal static class Util |
|||
{ |
|||
public static Tuple<TU[], TV[]> UnpackSinglePass<TU, TV>(this IEnumerable<Tuple<TU, TV>> samples) |
|||
{ |
|||
var u = new List<TU>(); |
|||
var v = new List<TV>(); |
|||
|
|||
foreach (var tuple in samples) |
|||
{ |
|||
u.Add(tuple.Item1); |
|||
v.Add(tuple.Item2); |
|||
} |
|||
|
|||
return new Tuple<TU[], TV[]>(u.ToArray(), v.ToArray()); |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,97 @@ |
|||
// <copyright file="LocallyWeightedRegression.cs" company="Math.NET">
|
|||
// Math.NET Numerics, part of the Math.NET Project
|
|||
// http://numerics.mathdotnet.com
|
|||
// http://github.com/mathnet/mathnet-numerics
|
|||
// http://mathnetnumerics.codeplex.com
|
|||
//
|
|||
// Copyright (c) 2009-2013 Math.NET
|
|||
//
|
|||
// Permission is hereby granted, free of charge, to any person
|
|||
// obtaining a copy of this software and associated documentation
|
|||
// files (the "Software"), to deal in the Software without
|
|||
// restriction, including without limitation the rights to use,
|
|||
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|||
// copies of the Software, and to permit persons to whom the
|
|||
// Software is furnished to do so, subject to the following
|
|||
// conditions:
|
|||
//
|
|||
// The above copyright notice and this permission notice shall be
|
|||
// included in all copies or substantial portions of the Software.
|
|||
//
|
|||
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|||
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
|||
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
|||
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
|||
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
|||
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|||
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
|||
// OTHER DEALINGS IN THE SOFTWARE.
|
|||
// </copyright>
|
|||
|
|||
using System; |
|||
using System.Collections.Generic; |
|||
using MathNet.Numerics.LinearAlgebra; |
|||
using MathNet.Numerics.LinearAlgebra.Storage; |
|||
|
|||
namespace MathNet.Numerics.LinearRegression |
|||
{ |
|||
public static class WeightedRegression |
|||
{ |
|||
/// <summary>
|
|||
/// Weighted Linear Regression using normal equations.
|
|||
/// </summary>
|
|||
public static Vector<T> Weighted<T>(Matrix<T> x, Vector<T> y, Matrix<T> w) where T : struct, IEquatable<T>, IFormattable |
|||
{ |
|||
return x.TransposeThisAndMultiply(w * x).Cholesky().Solve(x.Transpose() * (w * y)); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Weighted Linear Regression using normal equations.
|
|||
/// </summary>
|
|||
/// <param name="intercept">True if an intercept should be added as first artificial perdictor value. Default = false.</param>
|
|||
public static T[] Weighted<T>(T[][] x, T[] y, T[] w, bool intercept = false) where T : struct, IEquatable<T>, IFormattable |
|||
{ |
|||
var predictor = Matrix<T>.Build.DenseMatrixOfRowArrays(x); |
|||
if (intercept) |
|||
{ |
|||
predictor = predictor.InsertColumn(0, Vector<T>.Build.DenseVector(predictor.RowCount, Vector<T>.One)); |
|||
} |
|||
var response = Matrix<T>.Build.DenseVector(y); |
|||
var weights = Matrix<T>.Build.DiagonalMatrix(new DiagonalMatrixStorage<T>(predictor.RowCount, predictor.RowCount, w)); |
|||
return predictor.TransposeThisAndMultiply(weights * predictor).Cholesky().Solve(predictor.Transpose() * (weights * response)).ToArray(); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Weighted Linear Regression using normal equations.
|
|||
/// </summary>
|
|||
/// <param name="intercept">True if an intercept should be added as first artificial perdictor value. Default = false.</param>
|
|||
public static T[] Weighted<T>(IEnumerable<Tuple<T[], T>> samples, T[] w, bool intercept = false) where T : struct, IEquatable<T>, IFormattable |
|||
{ |
|||
var xy = samples.UnpackSinglePass(); |
|||
return Weighted(xy.Item1, xy.Item2, w, intercept); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Locally-Weighted Linear Regression using normal equations.
|
|||
/// </summary>
|
|||
public static Vector<T> Local<T>(Matrix<T> x, Vector<T> y, Vector<T> t, Func<Vector<T>, Vector<T>, T> kernel) where T : struct, IEquatable<T>, IFormattable |
|||
{ |
|||
// TODO: Kernel definition is a bit weird as it includes computing the difference norm
|
|||
// We can make this more common once we change the norm to always be of type double around LA.
|
|||
|
|||
var w = Matrix<T>.Build.DenseMatrix(x.RowCount, x.RowCount); |
|||
for (int i = 0; i < x.RowCount; i++) |
|||
{ |
|||
w.At(i, i, kernel(t, x.Row(i))); |
|||
} |
|||
return Weighted(x, y, w); |
|||
} |
|||
|
|||
public static Func<Vector<double>, Vector<double>, double> GaussianKernel(double radius) |
|||
{ |
|||
// TODO: see above...
|
|||
var d = -2.0*radius*radius; |
|||
return (t, x) => Math.Exp(Distance.SSD(x, t)/d); |
|||
} |
|||
} |
|||
} |
|||
Loading…
Reference in new issue