Using Least Squares to Fit a Plane to Noisy Data

I have generated a data set by taking the equation for a plane in the form z=Ax+By+C and adding noise to z with a random number generator, z=Ax+By+C+randomnoise(x,y). I did this as x & y take on integer values from 0 to 9. It is stored as a 10x10 matrix called Zplanedata. The goal is to find the best estimate for A, B and C using the least squares method.

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

First we read in the data file. It is in the form of a 10x10 matrix. The [i,j] entry is the z value at x=i-1, y=j-1.

>	read `Zplanedata1.m`;

Then I display the Matrix that was just read. This step is not necessary and would be impractical for very large data sets.

>	Zplanedata1;

Having read in the data matrix we store it in a big one dimensional array. First we make an empty array, then we use two nested do loops to transfer the data from Zplanedata1 to Zvector.

>	Zvector := array(1..100):

>	for i from 1 to 10 do   for j from 1 to 10 do      Zvector[(i-1)*10+j]:=Zplanedata1[i,j]; od; od;

Next we convert the data type of Zvector from array type to Vector type. This is a awkward step I wish I could find a way to avoid. The transpose is used to make Zvector a column vector instead of a row vector.

>	Zvector:=Transpose(convert(Zvector,Vector)):

Next we make the analog of the Vandermonde matrix. The first two columns track the grid locations for x and y. The third column is always 1 since the coefficient of C is always 1. (There is probably a better way to fill in the values for A than the method I used below.)

>	A:=array(1..100,1..3):

>	for j from 1 to 100 do A[j,1]:=(j-modp(j,10))/10+1;         A[j,2]:=modp(j,10); if(A[j,2]=0) then A[j,2]:=10 fi; A[j,3]:=1; if (A[j,2]=10)then A[j,1]:=A[j-1,1] fi; od:

>	A:=convert(A,Matrix):

We would like to solve Ax=Zvector, but unless every data point lies on the exact same plane, Ax=z vector will not have a solution. Instead we get a least squares estimate by solving A^T * A* x = A^T * Zvector for the 3x1 column vector x=[A,B,C]^T , although I called x the coeffvector below.

>	coeffvector:=LinearSolve(Transpose(A).A,Transpose(A).Zvector);

This is equivalent to

>	LeastSquares(A,Zvector);

These three numbers are the estimates for A, B and C, respectively.

We plot this plane, z=Ax+By+C, and the raw data together. First we create the plane plot and store it as plane.

>	plane:=plot3d((coeffvector[1])*x + (coeffvector[2])*y+coeffvector[3],x=0..9, y=0..9,color=red):

Before plotting the data we convert it to array format. (I'd like to find a way to avoid this.)

>	zdata:=convert(Zplanedata1,array):

Now we create the plot and store it.

>	noisydata:=plot3d(zdata[x+1,y+1],x=0..9,y=0..9,color=green):

Now we display the two plots together.

>	display(plane,noisydata,axes=normal);

By the way, the original plane was z=2x-3y+2.

Homework Problem: Read in Zplanedata2.m from http://galileo.math.siu.edu/221/Maple/Zplanedata2.m.

It was generated by adding noise to a plane given by z = Ax+By + C.
Use the methods above to estimate the equation used. Graph the data set and your plane. Clean up your worksheet, add an explanation of the two programs that were used to find create Zvector and the Vandermonde-like matrix A. Print and turn in. Be sure your name is printed on the page.