% OrthogPolyRegr: Orthogonal polynomial regression.
%
%     Usage: [b,R2a,F,pF,ssr,sse,ssto,dfr,dfe,dfto,msr,mse] = ...
%                 orthogpolyregr(x,y,degree,{doplot})
%
%       x,y =     vectors of independent and dependent variables.
%       degree =  degree of polynomial (=k).
%       doplot =  optional boolean flag indicating, if true, that a plot of data
%                   and fitted polynomial is to be produced [default=0].
%       -------------------------------------------------------------------------------
%       b =       [k+1 x 1]  vector of regression coefficients.
%       R2a =     matching vector of cumulative adjusted coefficients of determination.
%       F =       [k x 1] vector of incremental, independent F statistics.
%       pF =      corresponding vector of corresponding probabilities.
%       ssr =     [k x 1] vector of incremental, independent regression sums-of-squares.
%       sse =     error sum of squares for full model.
%       ssto =    total sum of squares.
%       dfr =     [k x 1] vector of regression df.
%       dfe =     error df.
%       dfto =    total df.
%       msr =     [k x 1] vector of regression mean squares.
%       mse =     error mean square.
%

% Dutka, AF and FJ Ewens. 1971. A method for improving the accuracy of polynomial
%   regression analysis. J Quality Technology 3:149-155.

% RE Strauss, 2/14/02

function [b,R2a,F,pF,ssr,sse,ssto,dfr,dfe,dfto,msr,mse] = orthogpolyregr(x,y,degree,doplot)
  if (nargin < 4) doplot = []; end;
  
  if (isempty(doplot))
    doplot = 0;
  end;

  x = x(:);
  y = y(:);
  n = length(x);
  if (n~=length(y))
    error('  ORTHOGPOLYREGR: data vectors are incompatible');
  end;
  
  [b,P,alpha,beta] = orthogpolyregrf(x,y,degree);
  
  ssr = zeros(degree+1,1);
  sse = zeros(degree+1,1);
  R2a = zeros(degree+1,1);
  pR2a = zeros(degree+1,1); 
  Fiof = zeros(degree+1,1);
  pFiof = zeros(degree+1,1);
  
  for i = 1:degree+1                      % Sums of squares
    [e,sse(i),ssr(i),ssto,R2a(i),pR2a(i),pred] = orthogpolyregre(y,b(1:i),P(:,1:i));
  end;
  sse = sse(end);
  ssr = ssr(2:end)-ssr(1:end-1);
  
  dfr = ones(degree,1);
  dfto = n-1;
  dfe = dfto - sum(dfr);
  
  mse = sse/dfe;
  msr = ssr./dfr;
  F = msr/mse;
  pF = 1-fcdf(F,dfr,dfe);
  
  if (doplot)
    plot(x,y,'ko');
    [xp,yp] = orthogpolyregrp(x,b,alpha,beta);
    hold on;
    plot(xp,yp,'k');
    hold off;
    putbnd(xp,yp);
  end;
  
  return;