%--------------------------------------------------------------------------
%function [h, outlier, outlier_index] = experiment_least_squares(A, b)
%--------------------------------------------------------------------------
%
% Matlab function to compute a least squares fit for data. This
% solves a least squares problem where the m x n matrix A has the
% m observations of n parameters as its rows and columns, resp.
% The vector b needs to be m x 1, consisting of the observed quantity.
% The residuals (amount each data point is off from the least squares
% fit) are plotted, once with the values and once with their absolute
% values.
%
% On return: 
%       h             = best fit parameters
%       outlier       = observation value that is worst fit
%       outlier_index = the row number ( = observation number) that gave
%                       the worst fit value.
%
%    Randall Bramley
%    Started: Wed 26 Mar 2008
%    Turned into a function: 26 Jan 2009
%    Modified: 11 Feb 2009, with comments.
%
%--------------------------------------------------------------
function [h, outlier, outlier_index] = experiment_least_squares(A, b)

    [number_of_observations, number_of_parameters] = size(A);
    
    %--------------------------------------------------
    % Solve the least squares problem min ||A*h - b||_2
    %--------------------------------------------------
    h = A\b;
    
    %--------------------------------------------------------------
    % Keep the print out of the answer for n = 3, until debugging
    % is finished.
    %--------------------------------------------------------------
    if (number_of_parameters == 3)
        s = sprintf('Best least squares fit for data is A = %g, B = %g , C = %g', ...
                  h(1), h(2), h(3));
    disp(s);
    end
    
    residual = b - A*h;
    display(['Residual vector is ', num2str(h')]);
    
    figure;
    plot(1:number_of_observations, residual, 'r-',...
         1:number_of_observations, residual, 'b*');
    title('Residual Vector');
    xlabel('Residual Vector Component');
    ylabel('Residual Vector Values');
    
    figure;
    plot(1:number_of_observations, abs(residual), 'r-',...
         1:number_of_observations, abs(residual), 'b*');
    title('Absolute Values of Residual Vector');
    xlabel('Residual Vector Component');
    ylabel('Residual Vector Values');

    unstack;
    
    [outlier, outlier_index] = max(abs(residual));

return