
clear all
close all force
  disp('Anacapa Data Bootstrap CIs')
  lw =2.5; 
  set(0, 'DefaultAxesFontSize', 16);
  fs = 16;
  msize = 10;
   s = RandStream('mt19937ar','Seed',0);
    RandStream.setGlobalStream(s);
    
    
%jabe = @(sample) mean(sample)/6 * (skewness(sample)^2 + (kurtosis(sample)-3)^2/4 );

muscaria=[10 11 12  9 10 11 13 12 10 11 ...
          11 13  9 10  9 10  8 12 10 11 ...
           9 10  7 11  8  9 11 11 10 12 ...
          10  8  7 11 12 10  9 10 11 10 ...
           8 10 10  8  9 10 13  9 12  9 ...
           9];


thetahat = mean(muscaria)   %9.5588
  
 B=20000;
 thetabs=[];
 for i = 1:B
   muscariab = bootsample(muscaria(:));
    thetabs=[thetabs mean(muscariab)];
 end

%Quantile Bootstrap CI
[quantile(thetabs,0.025)     quantile( thetabs,0.975)]   % 0.1685   42.7039
% Centered Bootstrap CI
[max(0,2*thetahat - quantile(thetabs,0.975) )    2*thetahat - quantile( thetabs,0.025)]
%    0   18.9491

%Bootstrap within bootstrap t-type CI
B1=500; B2=1000;
tbi=[]; means=[];
for i=1:B1 
    muscariab = bootsample(muscaria(:));
     samplej=[];
    for j=1:B2 %bootstrap-within-bootstrap
       muscariabb = bootsample(muscariab(:));
       samplej=[samplej mean(muscariabb)];
    end
     meanbar=mean(samplej); meanstd=std(samplej);
    tbi = [tbi (mean(muscariab) - meanbar)/meanstd];
    means=[means mean(muscariab)];
end
 
[mean(means) + quantile(tbi, 0.025)*std(means) ...
    mean(means) + quantile(tbi, 0.975)*std(means)]
    
 %1.4188   11.9744
