% ----------------------------------------------------------------------
% Linear regression of steam pressure with "polyfit"
% Programming Note: call "polyfit" and "polyval"
% Instructor: Nam Sun Wang
% ----------------------------------------------------------------------

% Start fresh ----------------------------------------------------------
      clear all

% Program header -------------------------------------------------------
      disp('Linear regression of steam pressure data with a polynomial')

% Read x-y data from a file --------------------------------------------
%     [ x y ] = read2col('steam.prn');
      load steam.dat;
      x = steam(:,1);
      y = steam(:,2);

      while 1
        iorder = input('Order of fit (quit=0, linear=2, quadratic=3, etc): ');
        if (iorder == 0) error('Quit by the user.'); end
        if (iorder < 0 ) error('Order must be positive.'); end

% Find the polynomial coefficients by calling the "polyfit" function ---
        a = polyfit(x, y, iorder-1);

% Print the coefficients (The coefficients are revered) ----------------
        disp('The coefficients are (from the lower order up): ')
        for i=1 : iorder
          fprintf('a(%2i) = %14.6e\n', i, a(iorder-i+1))
        end

% Generate the fitted polynomial function with "polyval" ---------------
        yfit = polyval(a, x);

% Print the residual error ---------------------------------------------
        sse = (yfit-y)' * (yfit-y);
        disp([ 'squared residual error = ', num2str(sse) ])

% Plot (every 5th data point) and compare results ----------------------
        plot(x(1:5:size(x,1)), y(1:5:size(y,1)), 'o', x, yfit, '-')
        xlabel('Temperature (F)')
        ylabel('Pressure (psia)')
        title('Steam Pressure')

      end