%-----------------------------------------------------
%
% Matlab script that uses linear least squares to find
% the heights of three hills from six observations.
% Those were specified as:
%  
% From ground level,
%   A is measured 1236 meters high
%   B is measured 1941 meters high
%   C is measured 2417 meters high
%  
% After climing  
%  
% From A, B is  711 meters higher 
% From A, C is 1177 meters higher 
% From B, C is  475 meters higher 
%
%-----------------------
%
% Started  : 11 Jan 2009
% Modified : 21 Jan 2009, for formatting and comments.
% Modified : 24 Feb 2009, to correct the last statement
%   of the last observation from A -> B to B -> C.
%
%-----------------------------------------------------

G = [ 1  0 0
      0  1 0
      0  0 1
      -1 1 0
      -1 0 1
      0 -1 1];

b = [1236; 1941; 2417; 711; 1177; 475];

%---------------------------------------------
% Get the number of parameters and experiments
%---------------------------------------------
[number_of_observations, number_of_parameters] = size(G);

%-----------------------------------------------------------
% Solve the least squares problem. Think of G\b as
% applying  inverse(G)*b. Since G is not a square matrix,
% it doesn't really have an inverse. Internally Matlab
% uses a "pseudo-inverse" approach about which we don't need
% to worry about the details. Do "help pinv" if you are
% curious about pseudo-inverses. Key thing here: once G and
% b are set up, it's a trivial one-liner to get h, the 
% vector of the three heights.
%-----------------------------------------------------------
h = G\b;

%----------------------------------------------------------------
% "residual" would be zero if all of the observations are exactly
% correct. So the size of it gives a measure of how accurate 
% they actually were.
%----------------------------------------------------------------
r = b - G*h;
disp('residual (error) vector = ');
disp(r);

s = ['Best least squares fit for data is height(A) = ', num2str(h(1)), ...
     ', height(B) = ', num2str(h(2)), ', height(C) = ', num2str(h(3))];
disp(s);
    
[outlier, outlier_index] = max(abs(r));
s = ['Worst outlier is observation number ', num2str(outlier_index)];
disp(s);