% OrthogPolyRegrF: Objective function for orthogonal polynomial regression
%
%     Usage: [b,ss,P,alpha,beta] = orthogpolyregrf(x,y,degree)
%
%         x =       column vector of unique independent-variable values.
%         y =    matching column vector of sums of dependent-variable values at 
%                     each level of x.
%         degree =  degree of polynomial.
%         ---------------------------------------------------------------------------
%         b =       column vector of regression coefficients.
%         P =       values of observations (rows) for each orthogonal polynomial term 
%                     (cols).
%         alpha =   alpha coefficients used to calculate polynomials.
%         beta =    beta coefficients used to calculate polynomials.

% 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,P,alpha,beta] = orthogpolyregrf(x,y,degree)
  k = length(x);
  
  P = zeros(k,degree+1);
  b = zeros(degree+1,1);
  ss = zeros(degree+1,1);
  alpha = zeros(degree+1,1);
  beta = zeros(degree+1,1);
  
  term = 0;                           % Constant term
  j = term+1;
  P(:,j) = ones(k,1);
  b(j) = sum(y.*P(:,j))./sum(P(:,j).^2);
  
  term = 1;                           % Linear term
  j = term+1;
  alpha(j) = sum(x.*P(:,j-1).^2)./sum(P(:,j-1).^2);   
  P(:,j) = (x-alpha(j)).*P(:,j-1);
  b(j) = sum(y.*P(:,j))./sum(P(:,j).^2);
  beta(j) = sum(x.*P(:,j-1).*P(:,j)) ./ sum(P(:,j-1).^2);

  for term = 2:degree                 % Remaining terms
    j = term+1;
    alpha(j) = sum(x.*P(:,j-1).^2)./sum(P(:,j-1).^2);
    P(:,j) = (x-alpha(j)).*P(:,j-1) - beta(j-1)*P(:,j-2);
    beta(j) = sum(x.*P(:,j-1).*P(:,j)) ./ sum(P(:,j-1).^2);
    b(j) = sum(y.*P(:,j))./sum(P(:,j).^2);
  end;

  return;