function [bE,rE,sE,PbE,CbE] = linregr(x,D,P,C,alpha)
%                 MULTIPLE LINEAR REGRESSION USING LEAST SQUARES
% for the linear model x = D*b + e, with random disturbances e: N(0,sigma^2*I)
% ----------------------------------------------------------------------------
% CALLING:
% 1) REGRESSION RESULT:   [bE,rE,sE,PbE,CbE] = linregr(x ,D,P,C,alpha)
% 2) REGRESSION OPERATOR: [R ,E0,1 ,PE0,CE0] = linregr([],D,P,C,alpha)
% INPUT:
% x  ... is an n x 1 vector of observations
%        x==[] => R,E0,PE0 and CE0 are matrices defining the regression operator
% D  ... is an n x p full-column rank matrix of regressors
% P  ... is an m x p matrix of regressors for prediction (optional:default=[])
% C  ... is an q x p matrix of parametric functions (optional:default=P)
% alpha ... is a percentual confidence risk (optional: default=[]=5%)
% OUTPUT:
% 1) REGRESSION RESULT:
% bE ... =[b,E] is an p x 2 matrix, where
%         b are parameter estimates
%         [b-E,b+E] are their confidence intervals of (100-alpha) %
% rE ... =[r,E] is an p x 2 matrix, where
%         r are residuals x-D*b
%         [r-E,r+E] are their confidence intervals of (100-alpha) %
%         (!!! E is not implemented in this version => rE is of size p x 1)
% sE ... =[s,E1,E2] is an 1 x 3 vector, where
%         s is the estimate of the standard deviation sigma
%         [s-E1,s+E2] is its confidence interval of (100-alpha) %
% PbE ...=[Pb,E] is an (m+1) x 2 matrix, where
%         Pb(1:m,:) are predicted values of x = P*b + e
%         Pb(m+1,:) is  predicted value  of y = sum(P)*b + sum(e)
%         [Pb-E,Pb+E] are their confidence intervals of (100-alpha) %
%         PbE = [] if P was not supplied on input (P=[])
% CbE ...=[Cb,E] is an q x 2 matrix, where
%         Cb are estimates of the parametric functions C*b
%         [Cb-E,Cb+E] are their confidence intervals of (100-alpha) %
%         CbE = [] if neither C nor P was supplied on input (C=P=[])
% 2) REGRESSION OPERATOR:
% R  ... p x n matrix defining the regression operator b = R*x
% E0 ... is an p x 1 matrix defining the confidence vector E from bE
%        by E=sigma.*E0 where sigma is an estimate of the standard deviation
%        of random disturbances, usually sigma=norm(D*b-x)/sqrt(n-p).
% sE ... sE = s = 1 ... unit standard deviation (scalar)
% PE0... is an (m+1) x 1 matrix determining the confidence matrix E from PbE
%        by E=sigma.*PE0 where sigma=norm(D*b-x)/sqrt(n-p) - see E0
% CE0... is an q x 1 matrix determining the confidence matrix E from CbE
%        by E=sigma.*CE0 where sigma=norm(D*b-x)/sqrt(n-p) - see E0

% (C) 1996, V. Vesely, Masaryk University of Brno (Czech Republic)
%
global ALPHA_linregr DF_linregr TP_linregr CHI_linregr

error(nargchk(2,5,nargin));
if nargin<3, P = []; end
if nargin<4, C = []; end
if nargin<5, alpha=[]; end
if isempty(alpha), alpha = 5; end

[n,p] = size(D);
alpha = 0.01*alpha;
alpha2 = alpha./2;
newalpha = isempty(ALPHA_linregr);                  % alpha-changed-flag
if ~newalpha, newalpha = (alpha~=ALPHA_linregr)|(n-p~=DF_linregr); end
if alpha<0 | alpha>100, error('Risk alpha is not in the range 0-100'); end
if ~isempty(x) & length(x) ~= n
   error('The length of x differs from the row size of D');
end
if p~=rank(D)
   error('The regressor matrix D must be of full-column rank');
end

[m,pp] = size(P);
if pp~=0 & p~=pp
   error('The matrices P and D must have equal column size');
end
[q,pp] = size(C);
if pp~=0 & p~=pp
   error('The matrices C and D must have equal column size');
end

if newalpha        % SAVE TIME AVOIDING REPEATED COMPUTATION OF PERCENTILES
                   % IN SUBSEQUENT CALLS
   tp    = tinv(1-alpha2,n-p);
           % (1-alpha/2)-percentile of the Student's T distribution
   chi1  = chi2inv(alpha2,n-p);
           % (alpha/2)-percentile of the chi-square distribution
   chi2  = chi2inv(1-alpha2,n-p);
           % (1-alpha/2)-percentile of the chi-square distribution
   ALPHA_linregr = alpha;
   DF_linregr    = n - p;
   TP_linregr = tp; CHI_linregr(1) = chi1; CHI_linregr(2) = chi2;
else
   tp = TP_linregr; chi1 = CHI_linregr(1); chi2 = CHI_linregr(2);
end

DDinv = inv(D'*D);
gamma = sqrt(diag(DDinv));
        % sigma*gamma(j) = standard deviation of the j-th parameter estimate
if ~isempty(x)
   x  = reshape(x,n,1);
   b  = DDinv*D'*x;        % parameter estimate
   r  = x - D*b;           % residuals
   rE = [r];               % r and half-widths of their confidence intervals
   Se = sum(r.^2);         % residual sum of squares
   s  = sqrt(Se./(n-p));   % estimate of sigma
   sE = [s,s-sqrt(Se./chi2),sqrt(Se./chi1)-s];
   bE = [b,tp*s*gamma];    % b and half-widths of their confidence intervals
else
   b  = [];
   bE = DDinv*D';
   s  = 1;                  % dummy estimate of sigma=1
   rE = tp*s*gamma;
   sE = s;
end
gamma = [];
if m>0 & nargout>=4        % If predictions required
   gamma = diag(P*DDinv*P');
   if m==1, Psum = P     ; gammas = 1 + gamma;
   else   , Psum = sum(P); gammas = m + Psum*DDinv*Psum';
   end
   P     = [P;Psum];       % Add sum of all as additional parametric function
   if isempty(b)
      PbE   = [tp*s*sqrt([1+gamma;gammas])];
   else
      PbE   = [P*b,tp*s*sqrt([1+gamma;gammas])];
   end
end
if nargout>=5              % If estimates for parametric functions required
   CbE = [];
   if q>0,
      if ~isempty(b), CbE = C*b; end
      gamma = diag(C*DDinv*C');
   elseif size(PbE,2)==2, CbE = PbE(1:m,1);
   end
   CbE = [CbE,tp*s*sqrt(gamma)];
end
