% Test the mean and median as measures of the central tendency of a Gaussian
% distribution.
function lorentzian_example(Ntrials)

t1 = cputime;  % record start time

N = 2.^(1:13)+1;
NN = numel(N)

for i=1:NN
   i
   Number_of_points = N(i)
   for trial=1:Ntrials
      data = cauchy_rnd(0,1,N(i),1);
      means(trial) = mean(data);
      medians(trial) = median(data);
   end
   
   means = sort(means);
   meanQ1(i) = means( round(0.25*Ntrials) );
   meanQ3(i) = means( round(0.75*Ntrials) );
   std_mean(i) = std(means);
   
   medians = sort(medians);
   medianQ1(i) = medians( round(0.25*Ntrials) );
   medianQ3(i) = medians( round(0.75*Ntrials) );
   std_median(i) = std(medians);
   
   t2 = cputime;
   time_elapsed_minutes = (t2-t1)/60
end


figure(1)
hold off
loglog(N, meanQ3-meanQ1, 'r-+')
hold on
title({'Convergence of sample mean and median', ...
   'for Cauchy distributed data, x_0=0, \gamma=1'})
xlabel('Number of data points')
ylabel(sprintf('Result of %i trials', Ntrials));
plot(N, medianQ3-medianQ1, 'b-+')
plot(N, std_mean, 'r--o')
plot(N, std_median, 'b--o')
legend('means Q3-Q1', 'medians Q3-Q1', 'means \sigma', 'medians \sigma')
legend('location','southwest')
orient tall
orient landscape
print('lorentzian_example_fig1.pdf', '-dpdf')

figure(2)
hold off
loglog(N, (meanQ3-meanQ1) ./ (medianQ3-medianQ1), 'k-+')
hold on
title({'Comparing statistical power of sample mean and median', ...
   'for Cauchy distributed data, x_0=0, \gamma=1'})
xlabel('Number of data points')
ylabel(sprintf('Result of %i trials', Ntrials));
plot(N, std_mean ./ std_median, 'k--o')
legend('(means Q3-Q1) / (medians Q3-Q1)','(means \sigma) / (medians \sigma)')
legend('location','northwest')
orient tall
orient landscape
print('lorentzian_example_fig2.pdf', '-dpdf')
