
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) length(sample)/6 * (skewness(sample)^2 + (kurtosis(sample)-3)^2/4 );

anacapa=[452   184   115   315   139   177   214   356   166   246   177   289   175 ... 
   296   205   324   260   188   208   109   204    89    320   256   138   198  ...
   191   193   316   122   305   203   396   250   230   214    46   256   204 ...
   150   218   261   143   229   173   132   175   236   220   212   119   144 ...
   147   171   216   232   216   164   185   216   199   236  237   206  87 ];



thetahat = jabe(anacapa)   %9.5588
% Unbiased 
% a=6*(n-2)/((n+1)*(n+3)); 
% b=3*(n-1)/(n+1); c= 24*n*(n-2)*(n-3)/((n+1)^2 *(n+3)*(n+5));
% chi2jbu =  skew^2/a + (kurt-b)^2/c
% 
 
 B=20000;
 thetabs=[];
 for i = 1:B
   anacapab = bootsample(anacapa(:));
    thetabs=[thetabs jabe(anacapab)];
 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=[]; jabes=[];
for i=1:B1 
    anacapab = bootsample(anacapa(:));
     samplej=[];
    for j=1:B2 %bootstrap-within-bootstrap
       anacapabb = bootsample(anacapab(:));
       samplej=[samplej jabe(anacapabb)];
    end
     jbbar=mean(samplej); jbstd=std(samplej);
    tbi = [tbi (jabe(anacapab) - jbbar)/jbstd];
    jabes=[jabes jabe(anacapab)];
end
 
[mean(jabes) + quantile(tbi, 0.025)*std(jabes) ...
    mean(jabes) + quantile(tbi, 0.975)*std(jabes)]
    
 %1.4188   11.9744
