%Hesterberg - Bootstrapping Example 18.7
%Bootstrap for Mean Difference - Verizon (ILEC) versus CLEC Service Times

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read Customer Service Times (hours) %
% N_Verizon=1664  N_ILEC = 23         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
fid=fopen([pwd,'\data\ca18_001.txt']);
fgetl(fid);  %skip header line
AoA=textscan(fid,'%d%s'); %read Numeric and Text columns from data file
stimes=AoA{1};   %extract service time data
providers=AoA{2}; %extract Service Provider labels
fclose(fid);

%segregate ILEC and CLEC service time data
% !!!!!!!!!! MODIFY CODE TO DUE THIS VIA THE PROVIDER TAG !!!!!!!!!!!!!
ilec=stimes(1:1664);
clec=stimes(1665:1687);

%seed the random number generator
rand('state',sum(100*clock));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use Resampling to generate Bootstrap Distribution %
% of Difference between Means                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
samplesize1=length(ilec);
samplesize2=length(clec);
n_resamples=10000;
for i=1:n_resamples
    %resample Verizon (ILEC) data
    rs1=randsample(ilec,samplesize1,true); %sample with replacement
    %resample Competing Carrier (CLEC) data
    rs2=randsample(clec,samplesize2,true);
    %
    remdiff(i)=mean(rs1)-mean(rs2); %calc/save resampled mean1-mean2
end
%
%compute and display summary statistics
remean=mean(remdiff);  %bootstrap distribution mean
restderr=std(remdiff); %bootstrap distribution standard error
remdiff=remdiff';
%percentile cutoffs
p1=prctile(remdiff,1);
p25=prctile(remdiff,2.5);
p5=prctile(remdiff,5);
p95=prctile(remdiff,95);
p975=prctile(remdiff,97.5);
p99=prctile(remdiff,99);
display('Bootstrap M1-M2 Distribution Summary statistics:');
display(['Mean:        ',num2str(remean)]);
display(['Std. Error:  ',num2str(restderr)]);
display(['1st  %ile:   ',num2str(p1)]);
display(['2.5  %ile:   ',num2str(p25)]);
display(['5th  %ile:   ',num2str(p5)]);
display(['95th %ile:   ',num2str(p95)]);
display(['97.5 %ile:   ',num2str(p975)]);
display(['99th %ile:   ',num2str(p99)]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot bootstrap distribution with parametric overlay %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear n; clear hours;
hours=[min(remdiff):1.0:max(remdiff)]; %data bins
n=histc(remdiff,hours);
figure(1); clf;
subplot(2,1,1);
set(gcf,'PaperPositionMode','auto');
set(gcf,'Position',[650 100 600 800]);
bar(hours,n,'k');
hold on;
axis([-25 10 0 max(n)+50]);
xlabel('Repair Time (hours)');
ylabel('Frequency');
title_str(1)={['Bootstrap Distribution of Difference between Means']};
title_str(2)={'Hesterberg Example 18.7'};
title_str(3)={['(Number of Resamples = ',num2str(n_resamples),')']};
title(title_str);
%annotate randomization distribution
f11str(1)={['Mean:        ',num2str(remean)]};
f11str(2)={['Std. Error:  ',num2str(restderr)]};
f11str(3)={['1st  %ile:   ',num2str(p1)]};
f11str(4)={['99th %ile:   ',num2str(p99)]};
text(0,max(n)/2,f11str);
%generate normal-quantile plot to assess normality
subplot(2,1,2);
qqplot(remdiff);
hold on;
ylabel('delta Repair Time (hours)');
title('');
hold off;

