function [a,y_fit,Sa,Sy_fit,Ra,R2,Vr,AIC,condNo] = mypolyfit(x,y,p,figNo,Sy,b)
% [a,y_fit,Sa,Sy_fit,Ra,R2,Vr,AIC,condNo] = mypolyfit(x,y,p,FigNo,Sy,b)
%
% fit a power-polynomial, y_fit(x;a) to data pairs (x,y) where
%
%  y_fit(x;a) = SUM_i=1^length(p)  a_i  x^p_i
%
% which minimizes the Chi-square error criterion, X2,
%
%  X2 = SUM_k=1 ^ length(x) [ ( y_fit(x_k;a) - y_k )^2 / Sy_k^2 ]
%
% where  Sy_k is the standard error of the k-th data point.  
%
%  INPUT     DESCRIPTION                                        DIM     DEFAULT
% ------     ------------------------------------------------  -----    -------
%   x        known vector of the independent variables,        (1xm)
%   y        measured vector of the dependent variables,       (1xm)
%   p        vector of the real power of x (x^p) of each term  (1xn)
%   figNo    figure number. Use 0 to suppress plotting                    0
%   Sy       measurement errors for each value of y  (1x1) or  (1xm)      1
%   b        regularization constant                                      0
%
% OUTPUT     DESCRIPTION                                         
% ------     ----------------------------------------------------
%   a        values of the polynomial coefficients
%   y_fit    poloynomial at values of x
%   Sa       standard errors of polynomial coefficients 
%   Sy_fit   standard errors of the curve-fit
%   Ra       parameter correlation matrix
%   R2       R-squared error criterion
%   Vr       unbiased estiamte of unweighted residuals
%   AIC      Akaike information criterion
%   condNo   condition number of the regularized system matrix 
%
% HP Gavin, 
% Fitting Models to Data: Generalized Linear Least Squares and Error Analysis
% https://people.duke.edu/~hpgavin/SystemID/linear-least-sqaures.pdf
%
% updated 2007-04-10, 2013-02-24

 if ( length(x) ~= length(y) )       % error checking
   disp(' the length of x must equal the length of y ');
   return
 end

 if nargin < 4,  FigNo = 0;  end     % default figure number
 if nargin < 5,  Sy = 1.0;   end     % default measurement error
 if nargin < 6,  b  = 0.0;   end     % default regularization parameter

 if any(Sy <= 0)                     % ensure that all vaues of Sy are positve
   Sy(find(Sy<=0)) = 1.0;
%  nargin = 4;                       % use Vr to compute invVy, Va, and Sy_fit
 end

 x = x(:);                           % make "x" a column-vector
 y = y(:);                           % make "y" a column-vector
 p = p(:);                           % make "p" a column-vector

 Nd = length(x);                     % number of data points
 Np = length(p);                     % number of parameters

 ISy = eye(Nd)./(Sy.^2);       % inverse of the measuremen error covariance, Vy

 X = zeros(Nd,Np);             % allocate memory for the model basis vectors

 for i=1:Np
   X(:,i) =  x.^p(i);          % the matrix of model basis vectors 
 end

 condNo = cond( X'*ISy*X + b*eye(Np) );              % condition number 

 a = inv( X'*ISy*X + b*eye(Np) ) * X'*ISy*y;         % least squares parameters

 y_fit = X*a;                                        % least squares fit

 Vr = sum((y-y_fit).^2) / (Nd-Np);   % unbiased variance of unweighted residuals

 if nargin < 5
   invVy = eye(Nd) / Vr;  % measurement error covariance computed from residuals
 else 
   invVy = ISy;           % measurement error covariance provided by user
 end

 Va = inv( X'*invVy*X + b*eye(Np) ); % estimated parameter covariance

 if b ~= 0                           % regularized least squares 
   Va = Va * (X'*ISy*Vy*ISy*X) * Va; % parameter covariance w/ regularization
 end

 Sa = sqrt(diag(Va));                % standard error of the parameters

 Ra = Va ./ (Sa * Sa');              % parameter cross-correlation matrix

 Sy_fit = sqrt(diag(X*Va*X'));       % standard error of the fit 

 R2  = 1 - sum( (y-y_fit).^2 ) / sum( (y-sum(y)/Nd).^2 ); % R-squared = correl

 AIC =  log(2*pi*Np*Vr) + (y_fit-y)'*invVy*(y_fit-y) + 2*Np; 

 disp('     p         a            +/-   da           (percent) ')
 disp('---------------------------------------------------------')
 for i=1:Np
   if rem(p,1) == 0
      fprintf('   a[%2d] =  %11.3e;    +/- %10.3e     (%6.2f %%)\n',  ...
                                p(i), a(i), Sa(i), 100*Sa(i)/abs(a(i)) ); 
   else
      fprintf(' %8.2f :  %11.3e     +/- %10.3e     (%6.2f %%)\n',  ...
                                p(i), a(i), Sa(i), 100*Sa(i)/abs(a(i)) ); 
   end
 end

if figNo % Plots ---------------------------------------
 
  formatPlot(12,2,5);  % font size, line width, marker size

  figure(figNo);        % plot the data, the model, and the confidence interval

    % confidence interval (CI) of the model (y_fit)
    yps95 =  y_fit + 1.96*Sy_fit;                    %  + 95 CI
    yms95 =  y_fit - 1.96*Sy_fit;                    %  - 95 CI
    yps99 =  y_fit + 2.58*Sy_fit;                    %  + 99 CI
    yms99 =  y_fit - 2.58*Sy_fit;                    %  - 99 CI
    xp  = [ x   ; x(end:-1:1)   ; x(1)   ];          % x coordinates for patch
    yp95  = [ yps95 ; yms95(end:-1:1) ; yps95(1) ];  % y coordinates for patch
    yp99  = [ yps99 ; yms99(end:-1:1) ; yps99(1) ];  % y coordinates for patch

    patchColor95 = [ 0.95, 0.95, 0.1 ];              % color for the 95% CI
    patchColor99 = [ 0.2,  0.95, 0.2 ];              % color for the 99% CI

    clf
    subplot(121)        % plot the data, curvefit, and CI's
    hold on
      hc99 = patch( xp, yp99, patchColor99);         % 99% CI
      hc95 = patch( xp, yp95, patchColor95);         % 95% CI
      hd   = plot( x, y,     'ob',  'LineWidth', 3); % data points
      hf   = plot( x, y_fit, '-k');                  % polynomial fit
    hold off
    set(hc95,'EdgeColor', patchColor95 );
    set(hc99,'EdgeColor', patchColor99 );
    xlabel('x') 
    legend([hd, hf , hc95, hc99 ], 'data', 'y_{fit}', '95% c.i.', '99% c.i.' );
    axis('tight')

    subplot(122)        % plot the data vs the model and show correlation 
    hold on 
      plot( y, y,'-k', 'LineWidth', 0.5 )             % 1-to-1 line
      plot( y_fit, y, 'ob', 'LineWidth', 3)           % data and model points
      tx = 0.90*min(y) + 0.10*max(y);
      ty = 0.02*min(y) + 0.98*max(y); 
      text(tx,1.0*ty, sprintf('condn # = %4.1f', condNo));
      text(tx,0.9*ty, sprintf('\\sigma_r = %5.3f', sqrt(Vr)));
      text(tx,0.8*ty, sprintf('AIC = %4.1f', AIC));
      text(tx,0.7*ty, sprintf('R^2 = %5.3f', R2));
      text(tx,0.6*ty, sprintf('n = %d', length(a) ));
    hold off
    xlabel('y_{fit}')
    ylabel('y')
    axis('tight')
% print('/home/hpgavin/Teaching/SystemID/CourseNotes/Figures/mypolyfit-1a.pdf','-dpdfcrop')

  figure(figNo+1)          % plot the histogram (distribution) of the residuals
    clf
    nBars = max(10,round(Nd/5)); 
    [fx,xx] = hist(y-y_fit, nBars, nBars/(max(x)-min(x)));
      bar(xx,fx)        % histogram of the residuals
      xlabel('residuals, r = y - y_{fit}')
      ylabel('empirical PDF, f_R(r)')
% print('/home/hpgavin/Teaching/SystemID/CourseNotes/Figures/mypolyfit-2a.pdf','-dpdfcrop')

  plotCDFci(y-y_fit, 95, figNo+2);  % plot the CDF with the CDF confidence intervals
  xlabel('residuals, r = y - y_{fit}')
  ylabel('CDF of residuals   F_R(r) = Prob[R <= r]')
% print('/home/hpgavin/Teaching/SystemID/CourseNotes/Figures/mypolyfit-3a.pdf','-dpdfcrop')

end % ------------------------------------------------------- Plots

% ------------------------------------------------------------------ mypolyfit
