close all force
clear all
%moon ilusion bootstapped
%
 disp('Moon Illusion Example in Bootstrap')
  lw = 2; 
  set(0, 'DefaultAxesFontSize', 16);
  fs = 16;
  msize = 10;
  rand('seed',0)

x =[1.73 1.06 2.03 1.40 0.95 1.13 1.41 1.73 1.63 1.56];
n=length(x);
mx =mean(x);
count=0;
t = (mx-1)/(std(x)/sqrt(n));
tstars =[];
B=100000
for b = 1:B
   xstar= bootsample(x');
   tstar= (mean(xstar) - mx)/(std(xstar)/sqrt(n));
   tstars =[tstars tstar];  
   count = count + (tstar > t);  %for p*-value
end
count/B  %p*value=0.0024
quantile(tstars, 1-0.05)  %1.8761
reject = t>quantile(tstars, 1-0.05)  %0 no; 1 yes; here 1.

figure(1)
histn(tstars, -5.5,0.25,5.5);
hold on
xx=-5.5:0.01:5.5;
plot(xx, tpdf(xx, 9),'r-','LineWidth',2)
plot(xx, normpdf(xx),'g-','LineWidth',2)
plot([t t],[0 0.45],'r-','LineWidth',2)
%add more
hold off
axis([-5.5 5.5 0 0.43])
%print -depsc 'C:\BESTAT\Boot\Booteps\moonillusionboot1.eps'

figure(2)
histn(tstars, -5.5,0.15,5.5);
hold on
xx=0:0.01:5.5;
plot(xx, tpdf(xx, 9),'r-','LineWidth',3)
plot(xx, normpdf(xx),'g-','LineWidth',3)
legend('bootstrap distribution','t','normal','Location','NorthEast')
plot([t t],[0 0.45],'r-','LineWidth',3)
%add more
hold off
axis([2.5 5.3 0 0.035])
%print -depsc 'C:\BESTAT\Boot\Booteps\moonillusionboot2.eps'

