
close all force
clear all
 disp(' Microdamage in Bones:  Bootstrap Approach')
  lw = 2; 
  set(0, 'DefaultAxesFontSize', 16);
  fs = 16;
  msize = 10;
  rand('seed',0)
  
%Microdamage in bones, Example 10.4,  page 365
x1=[0.790 0.944 0.958 1.011 0.714 0.256 0.406 ...
    0.135 0.316 1.264 1.410 1.160 0.179];
x2 =[1.137 0.601 1.029 1.264 1.183 1.856 1.899 ...
    0.486 0.813 1.327 1.325 2.012 1.026 1.130 0.605 0.870 0.820];

n1=length(x1); n2 = length(x2);
mx1 =mean(x1); mx2 = mean(x2);
 s1=std(x1); s2 = std(x2);
count=0;
 t = (mx1-mx2)/(sqrt(s1^2/n1 + s2^2/n2))
 ndf = (s1^2/n1 + s2^2/n2 )^2 /( (s1^2/n1)^2/(n1-1) + ...
     (s2^2/n2)^2 /(n2-1) )
%   Welch-Satterthwaite df = 26.4234
pval=tcdf(t, ndf)  %0.0096
tstars =[];
B=100000;
for b = 1:B
   x1star= bootsample(x1');
   x2star= bootsample(x2');
   mxs1 =mean(x1star); mxs2 = mean(x2star);
   ss1=std(x1star); ss2 = std(x2star);
   tstar= (mxs1 - mxs2 - (mx1 - mx2))/sqrt(ss1^2/n1 + ss2^2/n2);
   tstars =[tstars tstar];  
   count = count + (tstar < t);  %for p*-value
end
count/B  %p*value 0.0083
quantile(tstars,  0.05); %-1.6578 
     %compare to   tinv(0.05, 26.4234) = -1.7046
reject = t < quantile(tstars,  0.05)  %0 no; 1 yes; here = 1

figure(1)
histn(tstars, -4.5,0.25,4.5);
hold on
xx= -5:0.01:5;
plot(xx, tpdf(xx, ndf),'r-','Linewidth',2)
plot([t t],[0 0.45],'r-','LineWidth',2)
xlabel('t^*')
text(-2.4, 0.42, ' t', 'Color', 'k','FontSize',18);
axis([-4.5 4.5 0 0.45])
hold off
%print -depsc 'C:\STAT\Boot\Booteps\microdamageboot1.eps'

figure(2)
histn(tstars, -4.5,0.15,4.5);
hold on
xx= -4.5:0.01:4.5;
plot(xx, tpdf(xx, ndf),'r-','Linewidth',3)
plot(xx, normpdf(xx),'g-','Linewidth',3)
legend('bootstrap distribution','t','normal','Location','NorthWest')
plot([t t],[0 0.6],'r-','LineWidth',3)
xlabel('t^*')
hold off
axis([-4.5 -2.2 0 0.04])
%print -depsc 'C:\STAT\Boot\Booteps\microdamageboot2.eps'

