I used the MathNet.Iridium release because it is compatible with .NET 3.5 and VS2008. The method is based on the Vandermonde matrix. Then I created a class to hold my polynomial regression
using MathNet.Numerics.LinearAlgebra;
public class PolynomialRegression
{
Vector x_data, y_data, coef;
int order;
public PolynomialRegression(Vector x_data, Vector y_data, int order)
{
if (x_data.Length != y_data.Length)
{
throw new IndexOutOfRangeException();
}
this.x_data = x_data;
this.y_data = y_data;
this.order = order;
int N = x_data.Length;
Matrix A = new Matrix(N, order + 1);
for (int i = 0; i < N; i++)
{
A.SetRowVector( VandermondeRow(x_data[i]) , i);
}
// Least Squares of |y=A(x)*c|
// tr(A)*y = tr(A)*A*c
// inv(tr(A)*A)*tr(A)*y = c
Matrix At = Matrix.Transpose(A);
Matrix y2 = new Matrix(y_data, N);
coef = (At * A).Solve(At * y2).GetColumnVector(0);
}
Vector VandermondeRow(double x)
{
double[] row = new double[order + 1];
for (int i = 0; i <= order; i++)
{
row[i] = Math.Pow(x, i);
}
return new Vector(row);
}
public double Fit(double x)
{
return Vector.ScalarProduct( VandermondeRow(x) , coef);
}
public int Order { get { return order; } }
public Vector Coefficients { get { return coef; } }
public Vector XData { get { return x_data; } }
public Vector YData { get { return y_data; } }
}
which then I use it like this:
using MathNet.Numerics.LinearAlgebra;
class Program
{
static void Main(string[] args)
{
Vector x_data = new Vector(new double[] { 0, 1, 2, 3, 4 });
Vector y_data = new Vector(new double[] { 1.0, 1.4, 1.6, 1.3, 0.9 });
var poly = new PolynomialRegression(x_data, y_data, 2);
Console.WriteLine("{0,6}{1,9}", "x", "y");
for (int i = 0; i < 10; i++)
{
double x = (i * 0.5);
double y = poly.Fit(x);
Console.WriteLine("{0,6:F2}{1,9:F4}", x, y);
}
}
}
Calculated coefficients of [1,0.57,-0.15]
with the output:
x y
0.00 1.0000
0.50 1.2475
1.00 1.4200
1.50 1.5175
2.00 1.5400
2.50 1.4875
3.00 1.3600
3.50 1.1575
4.00 0.8800
4.50 0.5275
Which matches the quadratic results from Wolfram Alpha.
Edit 1
To get to the fit you want try the following initialization for x_data
and y_data
:
Matrix points = new Matrix( new double[,] { { 1, 82.96 },
{ 2, 86.23 }, { 3, 87.09 }, { 4, 84.28 },
{ 5, 83.69 }, { 6, 89.18 }, { 7, 85.71 },
{ 8, 85.05 }, { 9, 85.58 }, { 10, 86.95 },
{ 11, 87.95 }, { 12, 89.44 }, { 13, 93.47 } } );
Vector x_data = points.GetColumnVector(0);
Vector y_data = points.GetColumnVector(1);
which produces the following coefficients (from lowest power to highest)
Coef=[85.892,-0.5542,0.074990]
x y
0.00 85.8920
1.00 85.4127
2.00 85.0835
3.00 84.9043
4.00 84.8750
5.00 84.9957
6.00 85.2664
7.00 85.6871
8.00 86.2577
9.00 86.9783
10.00 87.8490
11.00 88.8695
12.00 90.0401
13.00 91.3607
14.00 92.8312