Least Squares Polynomial Fitting for Stored Data

In this worksheet we read in data from two files that I created earlier. Each was generated by adding random noise to a polynomial. The goal is to find the polynomial.

Down load the following files from my webpage:

http://galileo.math.siu.edu/Courses/221/Maple/datasetx.m

http://galileo.math.siu.edu/Courses/221/Maple/dataset3.m

Read in each data set using the read command. (You will have too figure out where they are stored on your computer.)

Each data set gives xdata, ydata, Ndata.

Ndata is the number of data points stored.

xdata is an array of length Ndata for the x-coordinates.

ydata is an array of length Ndata for the y-coordinates.

Example 1. We read in datasetx.m, plot it, find a least squares fit for a suitable polynomial, and finally plot the data with the curve.

>	restart:with(LinearAlgebra):with(plots):

Warning, the name changecoords has been redefined

>	read `datasetx.m`;

>

Ndata;

>	plot([[convert(xdata,list)[n],convert(ydata,list)[n]] $n=1..Ndata],style=POINT);

After plotting the data decide how many terms to use for a polynomial least sqaures fit. Then create the Vandermonde matrix and apply linsove. It looks like a cubic too me. So, we set Nterms to 4.

>	Nterms:=4:

>	V:=VandermondeMatrix(xdata,Ndata,Nterms):

>	a:=LinearSolve(Transpose(V).V,Transpose(V).ydata);

>	CurvePlot:=plot(sum('a[k]*x^(k-1)','k'=1..Nterms), x=xdata[1]..xdata[Ndata],color=black,thickness=3): PointPlot:=plot([[convert(xdata,list)[n],convert(ydata,list)[n]] $n=1..Ndata],style=POINT): display([PointPlot,CurvePlot]);

Looks pretty good! Lastly, we compute a messure of how good the fit it. It is just the square root of the sum of the squares of the differences between the curve and the ydata points.

>	sqrt(  sum( (sum(convert(a,list)[k]*convert(xdata,list)[n]^(k-1),k=1..Nterms)-convert(ydata,list)[n])^2,n=1..Ndata)  );

In fact this 30 point data set was from p(x)= -x^3+2x over [-2,2] with random noise between -1 and 1 added on.

Problem 1. Read in dataset3.m and do a least squares fit for a suitable polynomial. Plot the results.

>	restart:with(LinearAlgebra):with(plots):