MATLAB: MATLAB coding problem with least squares surface fit

least squares surface fitting

[EDIT: Thu May 19 17:56:20 UTC 2011 – Reformat – MKF]
Good (Proper Greeting of the Day All),
Again, I am very new to the MATLAB programming thing but I think I am getting better (yay!!) so any advice or help would be much appreciated. I am currently working on a program that will take data (x,y,z) and generate a surface from them using least squares approximation method. I am pretty sure this feature exists in one of MATLAB's amazing toolbars but I decided to take it upon myself and try to make my own. I want to plot the surface as well as give the error output. This is what I have so far but it says error:
??? Error using ==> mldivide
Matrix dimensions must agree.
Here is my coding:
function f=planefit(x,y,z)
%PLANEFIT(x, y, z) fits a least-squares polynomial through

%a set of data where x and y and z are the coordinates.

n = length(x); 
S1 = sum(x); S2 = sum(y); S3 = sum(z);
S4 = sum(x.*x); S5 = sum(x.*y); S6 = sum(y.*y);
S7 = sum(z.*x); S8 = sum(z.*y);
%Solving Matrix

A = [S4 S5 S1; S5 S6 S2; S1 S2 n];
B = [S7 S8 S3]; 
C = A\B;
%Surface Polynomial

for k=1:n
z = C(1)*x(k) + C(2)*y(k) + C(3); 
end
%Error Calculation

Error(k) = abs(z(k)-z);
Error = sum(Error.*Error);
%data visualization

[X Y}=meshgrid(x,y)
figure
surf(X,Y,z,'r')  %The Least Squares Surface

hold on
plot3(x, y,z 'o')  %The Original Data Points


grid on
EDIT After trying Sean de's suggestion.
Here is the new error I now have:
Error =
       0
??? Attempted to access z(2); index out of bounds because numel(z)=1.
Error in ==> planefit at 19 Error(k) = abs(z(k)-z)
Here is my current updated code:
function f=planefit(x,y,z)
%PLANEFIT(x, y, z) fits a least-squares polynomial through
%a set of data where x and y and z are the coordinates.
n = length(x);
S1 = sum(x); S2 = sum(y); S3 = sum(z);
S4 = sum(x.*x); S5 = sum(x.*y); S6 = sum(y.*y);
S7 = sum(z.*x); S8 = sum(z.*y);
%Solving Matrix
A = [S4 S5 S1; S5 S6 S2; S1 S2 n];
B = [S7; S8; S3];
C = A\B;
%Surface Polynomial
for k=1:n
    z = C(1)*x(k) + C(2)*y(k) + C(3);
    Error(k) = abs(z(k)-z)
end
%Error Calculation
Error = sum(Error.*Error)
%data visualization
[X Y]=meshgrid(x,y); figure
surf(X,Y,z,'r') %The Least Squares Surface
hold on
plot3(x, y,z,'o') %The Original Data Points
grid on
plot3(x, y,z,'o') %The Original Data Points
grid on

Best Answer

[EDIT: Thu May 19 18:44:56 UTC 2011 - Merged Comments - MKF]
Should B be a column vector? Either:
B = [S7;S8;S3];
or
B = [S7 S8 S3].';
EDIT
for k = 1:n 
   z(k) = stuff 
%more stuff 
end
right now you're overwriting z on each loop iteration.
Related Question