clear all
%close all
s = RandStream('mt19937ar','Seed',0);
RandStream.setGlobalStream(s);


 disp('Bootstrapping Regression: Body Fat Example')
 disp('Random Design Bootstrap')
 lw = 3; 
 set(0, 'DefaultAxesFontSize', 22);
 fs = 22;
 msize = 10;
load('C:\STAT\Reg\Regdat\fat.dat')
     casen = fat(:,1); %case number
broz = fat(:,2); %dependent variable Y
     siri = fat(:,3); %function of densi
     densi = fat(:,4); %invasive (makes people wet!)
age = fat(:,5); %below are predictors
weight = fat(:,6);   
height = fat(:,7);   
adiposi = fat(:,8); 
    % ffwei is functionally dependent on density//brozek and siri indices.
neck = fat(:,10); 
chest = fat(:,11);
abdomen  = fat(:,12);   
hip = fat(:,13); 
thigh = fat(:,14);
knee = fat(:,15);   
ankle = fat(:,16);
biceps = fat(:,17);
forearm = fat(:,18);     
wrist = fat(:,19);
vecones = ones(size(broz)); % necessary for intercept
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
 disp('=================================================================')
 disp('             p = 15, all 14 variables + intercept')
 disp('=================================================================')
% First to do standard analysis and get estimators of



X =[vecones age weight height adiposi  neck chest abdomen ...
    hip thigh knee ankle biceps forearm wrist];
y = broz;
% X is design matrix, n x p where n is the number of subjects and p is
% the number of parameters, or number of predictors + 1.
varnames = ['intercept=1 ' 'age=2 ' 'weight=3 ' ...
    'height=4 ' 'adiposi=5 ' 'neck=6 ' 'chest=7 ' 'abdomen=8 '...
    'hip=9' 'thigh=10 ' 'knee=11 ' 'ankle=12 ' 'biceps=13 ' 'forearm=14 ' 'wrist=15 '];
[n, p] = size(X)
% coefficients in multivariate linear regression
b = inv(X' * X) * X' * y 
% hat matrix for the predictions and residuals: yhat=H y.
H = X * inv(X' * X) * X';
yhat = H *y;
res = y - yhat;
%New person
Xh=[1  38  191  72  26  41   104  ...
   95  101.5   66  39   24   31   30  18.5];
%
  ypred = Xh  * b; %prediction of brozel index for 
  % a new person with covariates Xh.

 II = eye(n); %identity matrix
 JJ = ones(n); %matrix of ones, n x n
 SST = y' * (II - 1/n * JJ) * y
 SSE =  y' * (II - H) * y
 MSE = SSE/(n-p)
 %---------------------------------------
 Rsq = 1 - SSE/SST  %
 s = sqrt(MSE)
 
% BOOTSTRAP
all = [y X]; 
sbs=[];  Rsqbs=[];  yhatbs=[]; bbs=[];
B=50000  
for i=1:B
 newall = bootsample(all);
 yb = newall(:,1);
 newX = newall(:,2:end);
 bb = inv(newX' * newX) * newX' * yb ;
 bbs=[bbs bb(5)]; %hat beta_4 (coef for adiposity)
 Hb = newX * inv(newX' * newX) * newX';
  % Measurements for a new individual starting with 1 for 
  % the intercept
Xh=[1  38  191  72  26  41   104  ...
   95  101.5   66  39   24   31   30  18.5];
%predicted brozek index 19.5143
  Yh = Xh  * bb; 
 yhatbs=[yhatbs Yh]; 
 SSTb = yb' * (II - 1/n * JJ) * yb;
 SSEb =  yb' * (II - Hb) * yb;
 % of course SST = SSR + SSE-------------
 MSEb = SSEb/(n-p);
 %----------------------------------
 Rsqb = 1 - SSEb/SSTb;
 Rsqbs=[Rsqbs Rsqb];
 sb = sqrt(MSEb);
 sbs =[sbs sb];
end

close all
figure(1)
  histn(bbs, min(bbs), (max(bbs)-min(bbs))/100, max(bbs))
  hold on
  plot([b(5) b(5)],[0 1],'r-', 'LineWidth',2)
  xlabel('$\hat \beta_4$','fontsize', fs, 'Interpreter','LaTeX')
  axis([-3 5 0 1])
    print -depsc 'C:\STAT\NP\NPeps\fatregbootr1.eps'
%   subplot(222)
%   histn(yhatbs, min(yhatbs), (max(yhatbs)-min(yhatbs))/100, max(yhatbs))
%   hold on
%     plot([ypred ypred],[0 0.43],'r-', 'LineWidth',2)
%   xlabel('Prediction $\hat y_b$', 'Interpreter','LaTeX')
%   axis([15 24 0 0.43])
  figure(2)
  histn(Rsqbs, min(Rsqbs), (max(Rsqbs)-min(Rsqbs))/80, max(Rsqbs))
  hold on
    plot([Rsq Rsq],[0 23],'r-', 'LineWidth',2)
  xlabel('$R^2$', 'fontsize', fs,'Interpreter','LaTeX')
 axis([0.65 0.85 0 23])
    print -depsc 'C:\STAT\NP\NPeps\fatregbootr2.eps'
 figure(3)
  histn(sbs, min(sbs), (max(sbs)-min(sbs))/100, max(sbs))
  hold on
    plot([s s],[0 2.8],'r-', 'LineWidth',2)
  xlabel('$s$','fontsize', fs, 'Interpreter','LaTeX')
  axis([3.2 4.5 0 2.8])
   print -depsc 'C:\STAT\NP\NPeps\fatregbootr3.eps'

 mean(sbs)     %3.8649 
 mean(Rsqbs)   %0.7627
 mean(yhatbs)  %19.5640
 mean(bbs)    %0.3884
  
 median(sbs)    %3.8654
 median(Rsqbs)    %0.7635
 median(yhatbs)  %19.5637
 median(bbs)     %0.1977 

 var(sbs)      %0.0250
 var(Rsqbs)     %6.0129e-04
 var(yhatbs)   %1.0485
 var(bbs)     %0.5929

   


