Least Square Data Fitting

Given a set of data ((Subscript[x, 0]), (Subscript[y, 0])), ((Subscript[x, 1]), (Subscript[y, 1])) ..., ((Subscript[x, m]), (Subscript[y, m])), the least square fitting calls for finding a function that fits the data in the least square sense. To be more specific, let {(Subscript[g, 0])(x), (Subscript[g, 1])(x), ..., (Subscript[g, n])(x)} be a family of linearly independent functions, where we usually have m >> n (m is much larger than n), we want to choose the coefficients (Subscript[c, 0]), (Subscript[c, 1]), ..., (Subscript[c, n]) such that 
                    (∑_(i=0)^m)[(∑_(j=0)^n)(Subscript[c, j]) (Subscript[g, j])((Subscript[x, i])) - (Subscript[y, i])]^2  = (min_((Subscript[a, 0]),(Subscript[a, 1]),..., (Subscript[a, n]) ))(∑_(i=0)^m)[(∑_(j=0)^n)(Subscript[a, j]) (Subscript[g, j])((Subscript[x, i])) - (Subscript[y, i])]^2.
It is known that the coefficients (Subscript[c, 0]), (Subscript[c, 1]), ..., (Subscript[c, n]) satisfy the normal equation
                   (∑_(j=0)^n)(Subscript[c, j]) (∑_(i=0)^m) (Subscript[g, j])((Subscript[x, i])) (Subscript[g, k])((Subscript[x, i]))=(∑_(i=0)^m) (Subscript[g, k])((Subscript[x, i]))(Subscript[y, i]),    k = 0,1, ... , n.

Solving least square problem with Mathematica

Mathematica has a built in function called Fit. Here are some examples:

Examle 1:  A linear fit:

Clear[x,XY];
XY={{1, 3.7}, {1.5, 3.75},{2, 4.1}, {2.2, 4.6},{2.5, 4.3}, {3, 5}};
f[x_]=Fit[XY, {1,x}, x]

dots = ListPlot[XY,PlotStyle->{PointSize[0.02]},DisplayFunction-> Identity];
least= Plot[f[x],{x,1,3.2},DisplayFunction->Identity];
Show[dots,least,  DisplayFunction->$DisplayFunction, PlotLabel->txt,AxesLabel->{"x","y"}]

Examle 2:  By quadratic polynomial:

XY={{0.24, 0.23}, {0.65, -0.26}, {0.95, -1.10},{1.24, -0.45}, {1.73,0.27},  {2.01, 0.10}, {2.23, -0.19},{2.52, 0.24}, {2.77,0.56},  {2.99, 1.00}};
f[x_]=Fit[XY, {1,x,x^2}, x]

dots = ListPlot[XY,PlotStyle->{PointSize[0.02]},DisplayFunction-> Identity];
least= Plot[f[x],{x,0,3.2},DisplayFunction->Identity];
Show[dots,least,DisplayFunction->$DisplayFunction, PlotLabel->txt,AxesLabel->{"x","y"}]

By polynomial of degree 6:

f[x_]=Fit[XY, {1,x,x^2, x^3, x^4, x^5,x^6}, x]

dots = ListPlot[XY,PlotStyle->{PointSize[0.02]},DisplayFunction-> Identity];
least= Plot[f[x],{x,0.1,3.2},DisplayFunction->Identity];
Show[dots,least,DisplayFunction->$DisplayFunction, AxesLabel->{"x","y"}]

Solving least square problem with orthogonal polynomias

If we choose {(Subscript[g, 0])(x), (Subscript[g, 1])(x), ..., (Subscript[g, n])(x)} be a basis for polynomials of degree at most n, then the best choice for the least square fitting is to choose them as orthonormal polynomials with respect to the inner product 
                     < f, g> = (∑_(j=0)^m)f((Subscript[x, i])) g((Subscript[x, i])),  
that is, we choose (Subscript[g, 0])(x), (Subscript[g, 1])(x), ..., (Subscript[g, n])(x) such that  < (Subscript[g, i]), (Subscript[g, j])> = 0 if i≠j and  < (Subscript[g, i]), (Subscript[g, i])> = 1, with  (Subscript[g, i])(x) being a polynomial of degree i. Then, if  (Subscript[y, i])= f((Subscript[x, i])) for i=0,1,..., m, then the least square polynomial is given by 
                      p(x) = (∑_(j=0)^n)<f,(Subscript[g, j])>(Subscript[g, j])(x).
The orthonormal polynomials are generated recursively (Gram-Schmidt process).

Algorithm

BestFit[XY0_,p0_]:=Block[{XY},XY=XY0;p=p0; m=Length[XY];
Clear[a,b,ips,P];
ips[f_,g_]:=Sum[f g/.x->XY[[i,1]],{i,1,m}];a[n_]:=a[n]=ips[x P[n-1,x],P[n-1,x]]/ips[P[n-1,x],P[n-1,x]];
b[n_]:=b[n]=ips[x P[n-1,x],P[n-2,x]]/ips[P[n-2,x],P[n-2,x]];
P[0,x_]=1;P[1,x_]=x-a[1]; P[n_,x_]:=P[n,x]=(x-a[n])P[n-1,x]-b[n] P[n-2,x];
For[i=1,i≤m,i++,y[XY[[i,1]]]=XY[[i,2]]];
BestPoly[x_,p]=Sum[(ips[y[x],P[k,x]]/ips[P[k,x],P[k,x]]) P[k,x],{k,0,p}];];

The function of the algorithm is of the form BestFit[XY, p], where XY is the data XY =  {((Subscript[x, 0]), (Subscript[y, 0])), ((Subscript[x, 1]), (Subscript[y, 1])) ..., ((Subscript[x, m]), (Subscript[y, m]))} and p is the degree of least square polynomial. The output is : orthogonal polynomials P[n,x], and least square polynomial BestPoly[x,p].

Example

Examle 1:  Use polynomials of degree 6:

XY={{0.24, 0.23}, {0.65, -0.26}, {0.95, -1.10},{1.24, -0.45}, {1.73,0.27}, {2.01, 0.10}, {2.23,-0.19}, {2.52, 0.24}, {2.77,0.56}, {2.99, 1.00}};  
BestFit[XY,6] 
Simplify[BestPoly[x,6]]

Here are the list of orthogonal polynomials up to degree 6:

For[j=0, j<=6, j++,
Print["P",j,"(x)=" Simplify[P[j,x]]] ]

nodes=ListPlot[XY,PlotStyle->PointSize[.02], DisplayFunction->Identity]; 
fcn=Plot[ BestPoly[x,6],{x,0.2,3.2}, DisplayFunction->Identity];
Show[nodes,fcn,DisplayFunction-> $DisplayFunction]

Examle 2:  Use polynomials of degree 6 on the data set defined via f(x) =  1 + sin (10 x):

g[x_]:=1 + Sin[10 x]; pt[x_]:= 10+(x-1)/14; n=15;
XY= Table[{ pt[i], N[g[pt[i]],10]}, {i,1,n}];

BestFit[XY,6]; 
Simplify[BestPoly[x,6]]

nodes=ListPlot[XY,PlotStyle->PointSize[.02], DisplayFunction->Identity];  
best=Plot[ BestPoly[x,6],{x,10,11}, DisplayFunction->Identity];
Show[nodes,best,DisplayFunction-> $DisplayFunction]

Comparing the result with the Mathematica fitting:

f[x_]=Fit[XY, {1,x,x^2, x^3, x^4, x^5,x^6}, x]

dots = ListPlot[XY,PlotStyle->{PointSize[0.02]},DisplayFunction->Identity];
least= Plot[f[x],{x,10,11}, DisplayFunction->Identity];
Show[dots,least,DisplayFunction-> $DisplayFunction,  AxesLabel->{"x","y    "}]

The comparison shows that the method via orthogonal polynomial gives a much better result. The reason is that the normal equation with the basis {1, x,  (x^2), (x^3),  (x^4), (x^5), (x^6)} is severly ill-conditioned in this case.