%Coagulation times by bootstrap
%Here 24 animals (Box, Hunter, Hunter; Statistics for Experimenters,
%p. 166) are randomly allocated to 4 different diets,
%but the numbers allocated to different diets are not the same.
%The coagulation time for blood is measured for each animal.
%Are the diet-based differences significant?
clear all
close all force
  disp('Coagulation Times by Bootstrap')
  lw =2.5; 
  set(0, 'DefaultAxesFontSize', 16);
  fs = 16;
  msize = 10;
   s = RandStream('mt19937ar','Seed',0);
    RandStream.setGlobalStream(s);
    


 times = [62, 60, 63, 59, ...
                63, 67, 71, 64, 65, 66, ...
                68, 66, 71, 67, 68, 68, ...
                56, 62, 60, 61, 63, 64, 63, 59];
 diets = [1 1 1 1   2 2 2 2 2 2   3 3 3 3 3 3  4 4 4 4 4 4 4 4];

ni =[4 6 6 8]; N=24; k=4;
barxi = [mean(times(diets==1))  mean(times(diets==2))  ...
               mean(times(diets==3))  mean(times(diets==4)) ];
s2i = [var(times(diets==1))  var(times(diets==2))  ...
               var(times(diets==3))  var(times(diets==4)) ];
barx = sum(ni.*barxi)/N;
Fnum= 1/(k-1) * sum( ni .* (barxi - barx).^2 );
Fdenom = 1/(N-k) * sum( (ni - 1) .* s2i );
F=Fnum/Fdenom;

B=100000;
Fstar=[];
for i= 1:B
select =  [randsample(1:4, 4, 'true')  randsample(5:10, 6, 'true')  ...
    randsample(11:16, 6, 'true')  randsample(17:24, 8, 'true')];
timesstar = times(select);
barxistar = [mean(timesstar(diets==1))  mean(timesstar(diets==2))  ...
               mean(timesstar(diets==3))  mean(timesstar(diets==4)) ];
s2istar = [var(timesstar(diets==1))  var(timesstar(diets==2))  ...
               var(timesstar(diets==3))  var(timesstar(diets==4)) ];
barxstar = sum(ni.*barxistar)/N; 
Fnumstar= 1/(k-1) * (sum( ni .* (barxistar - barxi).^2)   - N*(barxstar-barx)^2  );
Fdenomstar = 1/(N-k) * sum( (ni - 1) .* s2istar );
Fstar=[Fstar Fnumstar/Fdenomstar];
end
histn(Fstar, 0, 0.2, 30);
hold on
plot(0:0.01:30, fpdf(0:0.01:30, k-1, N-k),'r-','Linewidth',2)
plot([F,F],[0, 0.9],'g-','LineWidth',2)
axis([0 20 0 1])
legend('bootstrap F*','pdf F_{3, 20}','F')
1- sum(F > Fstar)/B  % 
%print -depsc 'C:\STAT\Boot\Booteps\coagulationanovaboot.eps'


           
  [p,table,stats] = anova1(times, diets)          
