%Hesterberg - Bootstrapping Sample Mean Example 18.2

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read random sample of Verizon Customer Service Times (hours) %
% N=1664                                                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
load([pwd,'\data\eg18_002.txt']);    %fetch data from subdirectory
verizon = eg18_002; clear eg18_002;  %rename the dataset to 'verizon'
verizonmean=mean(verizon); %sample mean
verizonstderr=std(verizon)/sqrt(length(verizon)); %standard error

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute and plot Verizon data histogram %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
hours=[0:max(verizon)]; %define histogram categories
n=histc(verizon,hours); %count events per histogram category
figure(1); clf;
subplot(2,1,1);
set(gcf,'PaperPositionMode','auto');
set(gcf,'Position',[30 100 600 600]);
bar(hours,n,'k');
hold on;
xlabel('Repair Time (hours)');
ylabel('Frequency');
title('Verizon Customer Service Repair Times');
%annotate empirical data sample
f11str(1)={['  Mean: ',num2str(verizonmean)]};
f11str(2)={['  SE:   ',num2str(verizonstderr)]};
f11str(3)={['  N:    ',num2str(length(verizon))]};
text(100,max(n)/2,f11str);
hold off;
%generate normal-quantile plot to assess normality
subplot(2,1,2);
qqplot(verizon);
hold on;
ylabel('Repair Time (hours)');
title('');
hold off;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use Resampling to generate Bootstrap Distrimution of the Mean %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
samplesize=length(verizon);
n_resamples=5000; %number of resamples
for i=1:n_resamples
    rs=randsample(verizon,samplesize,true); %sample with replacement
    remeans(i)=mean(rs); %calc/save resampled mean
end
remean=mean(remeans);  %bootstrap distribution mean
restderr=std(remeans); %bootstrap distribution standard error
display('Bootstrapping Mean Summary:');
display('Original Data');
display(['  Mean:   ',num2str(verizonmean)]);
display(['  SE:     ',num2str(verizonstderr)]);
display('Bootstrap Distribution');
display(['  Mean:   ',num2str(remean)]);
display(['  Bias:   ',num2str(verizonmean-remean)]);
display(['  SE:     ',num2str(restderr)]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute parametric sampling distribution of the mean %
% for comparison to the bootstrap distribution         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sdm_hours=[verizonmean-(4*verizonstderr):0.01:verizonmean+(4*verizonstderr)];
for i=1:length(sdm_hours);
    sdm_z(i)=(sdm_hours(i)-verizonmean)/verizonstderr;
    sdm_pdf(i)=normpdf(sdm_z(i),0,1)/normpdf(0,0,1);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot bootstrap distribution with parametric overlay %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear n; clear hours;
hours=[6.0:0.05:10.0]; %data bins
n=histc(remeans,hours);
figure(2); clf;
subplot(3,1,1);
set(gcf,'PaperPositionMode','auto');
set(gcf,'Position',[650 100 600 800]);
bar(hours,n,'k');
hold on;
xlabel('Repair Time (hours)');
ylabel('Frequency');
titlestr(1)={'Bootstrap Distribution of Resampled Means'};
titlestr(2)={'Hesterberg Example 18.2'};
titlestr(3)={['Number of Resamples = ',num2str(n_resamples)]};
title(titlestr);
%overlay parametric Sampling Distribution of the Mean
%but first scale pdf to fit existing y-axis
sdm_pdf=sdm_pdf*max(n);
plot(sdm_hours,sdm_pdf,'r-');
legend('Bootstrap','SDM','Location','NorthEast');
%annotate bootstrap distribution
f21str(1)={['  Mean: ',num2str(remean)]};
f21str(2)={['  Bias: ',num2str(verizonmean-remean)]};
f21str(3)={['  SE:   ',num2str(restderr)]};
f21str(4)={[' 1%ile: ',num2str(prctile(remeans,1))]};
f21str(5)={['99%ile: ',num2str(prctile(remeans,99))]};
text(6.0,max(n)/2,f21str);
hold off;
%generate normal-quantile plot to assess normality
subplot(3,1,2);
qqplot(remeans);
hold on;
ylabel('Repair Time (hours)');
title('');
hold off;
%generate and plot cummulative empirical probability distribution
subplot(3,1,3);
[h,stats]=cdfplot(remeans);
hold on;
xlabel('Repair Time (hours)');
ylabel('Cummulative Probability');
title('');
hold off;
