function [jkResSs] = olsJackknife(X, y)
%olsJackknife Computes Jackknife residual sum of squares for a multiple regression.
%
% This program computes the Jackknife residual sum of squares for a multiple regression
% problem.  A jackknife parameter estimate is computed for a data set without i-th
% observation, i = 1, 2, ..., n, where n is the total number of observations.
% Estimated y value for the i-th observation is calculated.
% The residual for the jackknife estimate (= yhat - ytrue), then, is computed.
%
% USE:
% [jackknifeSS] = olsJackknife(X, y)
% 
% where X is the data matrix (may contain the intercept term),
% y is the dependent variable vector,
% jackknifeSS = the sum of jackknife residual sum of squares.
%
% It uses olsEstimates function.  

% Tomo Eguchi
% 28 January 2002

%X1 = normrnd(4,3,30,1);X2=normrnd(5,20,30,1);X3=normrnd(12,3,30,1);X4=normrnd(0,3,30,1);
%X = [ones(30, 1), X1, X2, X3, X4];
%trueA = [4; 3; 12; 9; 5];
%y = X*trueA + normrnd(0, 1, 30, 1);

[nrX, ncX] = size(X);             % figure out the size of the data matrix
jkResSs = 0;                      % initialize the output 

for i = 1:nrX,                    % for each row of the data matrix
   JKData = zeros(nrX-1, ncX);    % take out one row at a time.
   JKy = zeros(nrX-1, 1);         % create JKdata and JKy
   if i == 1,
      JKData = X(2:end,:);
      JKy = y(2:end);
   elseif i == nrX,
      JKData = X(1:(end-1),:);
      JKy = y(1:(end-1));
   else
      JKData = [X(1:i-1,:); X(i+1:end,:)];
      JKy = [y(1:i-1); y(i+1:end)];
   end
   
   [betaHat, aHat] = olsEstimates(JKData, JKy);   % compute the parameters for these.
   % error checking here:
%   [aHat] = regress(JKy, JKData);       % using the matlab regress command.
%   sum(aHat - betaHat2)                     % this should be very close to zero

   if (isnan(betaHat) == 1),                      % if NaN, flag.
      fprintf('Beta hat is NaN at i = %d. --- olsJackknife.m\n', i); 
   end
   
   yHat = X(i,:) * aHat;                          % compute yHat for the left-out datum.
   jkResSs = jkResSs + sum((yHat - y(i)).^2);      % compute the residual sum of squares
end