Hi M,
Here is a way to fit a line in 2d with the equation
a1*x1 + a2*x2 = 1
where (x1,x2) are (x,y). This representation gets rid of infinite slope problems.
The same code works in general in m dimensions to fit an m-1 dimensional plane
a1*x1 + a2*x2 + ... am*xm = 1.
For n points, let X be an nxm matrix where each row in X represents a point in m dimensions.
X = rand(20,2) % for example
Xcm = mean(X); % coordinates of center of mass
[~,s,v] = svd(X-Xcm,'econ');
n = v(:,end); % unit vector for the smallest singular value (they are sorted)
L = Xcm*n; % displacement of the best fit plane from the origin
if L < 0 % turn displacement into distance
L = -L;
n = -n;
end
a = n/L; % fitted plane, a1*x1 + a2*x2 + ... am*xm = 1.
drms = rms(X*n -L);
% L = 1/norm(a), so drms also equals (1/norm(a))*rms(X*a-1)
drms is the rms distance of the set of points to the plane.
For older versions of Matlab, in the svd line you would have to use
X-repmat(Xcm,size(X,1),1) in place of X-Xcm
