%Hesterberg - Bootstrapping Example 18.10
%Correlation of Baseball salaries and batting averages

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read original sample data   %
% Name|Salary|Average         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
fid=fopen([pwd,'\data\eg18_010.txt']);
fgetl(fid);  %skip header line
AoA=textscan(fid,'%s%f%f','delimiter','\t');
players=AoA{1};
salary=AoA{2};  %extract salary column
batting=AoA{3}; %extract batting average column
fclose(fid);

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use Resampling to generate Bootstrap Distribution %
% of the correlation coefficient                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
samplesize=length(salary);
indices=[1:samplesize]; %proxy indices for baseball players' data pair
n_resamples=10000;
for k=1:n_resamples
    %resample with replacement
    ri=randsample(indices,samplesize,true)';
    rs=salary(ri);
    rb=batting(ri);
    %compute/save correlation coefficient
    r=corr([rs rb]);
    recorr(k)=r(1,2);
end
%
%compute and display summary statistics
recorrmean=mean(recorr);  %bootstrap distribution mean
recorrstderr=std(recorr); %bootstrap distribution standard error
%
%percentile cutoffs
p1=prctile(recorr,1);
p25=prctile(recorr,2.5);
p5=prctile(recorr,5);
p95=prctile(recorr,95);
p975=prctile(recorr,97.5);
p99=prctile(recorr,99);
display('Correlation Bootstrap Summary statistics:');
display(['Mean Corr:   ',num2str(recorrmean)]);
display(['Std. Error:  ',num2str(recorrstderr)]);
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                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x=[min(recorr):0.05:max(recorr)]; %data bins
n=histc(recorr,x);
figure(1); clf;
subplot(3,1,1);
set(gcf,'PaperPositionMode','auto');
set(gcf,'Position',[650 100 600 800]);
bar(x,n,'k');
hold on;
axis([-0.6 0.6 0 max(n)+50]);
xlabel('Correlation Coefficient');
ylabel('Frequency');
title_str(1)={['Bootstrap Distribution of Correlation Coeficient']};
title_str(2)={['Hesterberg Example 18.10']};
title_str(3)={['(Number of Resamples = ',num2str(n_resamples),')']};
title(title_str);
%annotate bootstrap distribution
f21str(1)={['Mean:  ',num2str(recorrmean)]};
f21str(2)={['SE  : ',num2str(recorrstderr)]};
f21str(3)={[' 1st %ile:   ',num2str(p1)]};
f21str(4)={['99th %ile:   ',num2str(p99)]};
text(-0.45,max(n)/2,f21str);
%generate and display normal-quantile plot
subplot(3,1,2);
qqplot(recorr);
hold on;
ylabel('Correlation Coefficient');
title('');
hold off;
%generate and plot cummulative empirical probability distribution
subplot(3,1,3);
[h,stats]=cdfplot(recorr);
hold on;
xlabel('Correlation Coefficient');
ylabel('Cummulative Probability');
title('');
hold off;
