% Test the ability of the bootstrap to evaluate the usefulness of 
% mean vs. median for particular data sets drawn from the 
% Cauchy-Lorentz distribution.
function lorentzian_bootstrap(Npoints)

t1 = cputime;  % record start time

N_bootstrap = 10000
N_datasets = 100

for i=1:N_datasets
   data_orig = cauchy_rnd(0, 1, Npoints, 1);
   mean_orig(i) = mean(data_orig);
   median_orig(i) = median(data_orig);
   
   for trial=1:N_bootstrap
      data = data_orig( ceil( rand(Npoints,1) * Npoints ) );
      means(trial) = mean(data);
      medians(trial) = median(data);
   end

   means = sort(means);
   meanPlo(i) = means( round(0.01*N_bootstrap) );
   meanPhi(i) = means( round(0.99*N_bootstrap) );
   std_mean(i) = std(means);

   medians = sort(medians);
   medianPlo(i) = medians( round(0.01*N_bootstrap) );
   medianPhi(i) = medians( round(0.99*N_bootstrap) );
   std_median(i) = std(medians);
   
   printf('Dataset %i completed at t=%f seconds.\n', i, (cputime-t1) )
end




figure(1)
hold off
loglog( abs(mean_orig), (meanPhi-meanPlo)/2, 'r+' )
hold on
Npoints_string = sprintf('%i', Npoints);
N_datasets_string = sprintf('%i', N_datasets);
title(['Bootstrap of sample mean and median for ',N_datasets_string,' datasets'; ...
   'Npoints=',Npoints_string,', Cauchy distributed data, x_0=0, \gamma=1'])
xlabel('Absolute error of statistic for original dataset')
ylabel('98% confidence limit of statistic based on bootstrap');
plot( abs(median_orig), (medianPhi-medianPlo)/2, 'b+')
foo = sort(abs(mean_orig));
plot( foo, foo, 'k-' )
   % y=x line for reference
legend('means (Phi-Plo)/2', 'medians (Phi-Plo)/2', 'err = confidence limit')
legend('location','northwest')
orient tall
orient landscape
print(['lorentzian_bootstrap_Npts',Npoints_string,'_fig1.pdf'], '-dpdf')


figure(2)
hold off
loglog( (meanPhi-meanPlo)/2, (medianPhi-medianPlo)/2, 'k+' )
hold on
title(['Bootstrap of sample mean and median for ',N_datasets_string,' datasets'; ...
   'Npoints=',Npoints_string,', Cauchy distributed data, x_0=0, \gamma=1'])
xlabel('mean')
ylabel('median');
plot( abs(mean_orig), abs(median_orig), 'm+' )
foo = sort( (meanPhi-meanPlo)/2 );
plot( foo, foo, 'k--')
legend('98% confidence interval (half-width)', 'absolute error', 'median = mean')
legend('location','southeast')
orient tall
orient landscape
print(['lorentzian_bootstrap_Npts',Npoints_string,'_fig2.pdf'], '-dpdf')
