%
% epage 325
% 
% Written by:
% -- 
% John L. Weatherwax                2008-02-20
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

close all; drawnow; 
clc; 
clear; addpath( '../../Code/eda_data' ); addpath( '../../Code/eda_toolbox' ); addpath( '../Chapter1' ); 

%
% create four plots using each of the rules for histogram bin widths ... for both data sets: 
%
for pi=1:2, 

  if( pi==1 )
    load galaxy; x=EastWest; n=length(x);
  else
    load forearm; x=forearm; n=length(x);
  end
  
  % get the limits, bins, bin centers etc:
  x_lim_left = min(x)-1; 
  x_lim_rght = max(x)+1;

  t0   = x_lim_left;
  tm   = x_lim_rght; 
  rng  = tm-t0; 


  % 1) use Sturge's rule: 
  % 
  nbin = ceil( 1 + log2(n) ); 
  hS   = rng/nbin;               % <- compute the width depending on the number of bins
  h    = hS; 
  bins = t0:h:(nbin*h+t0);       % <- the bin edges ... 
  bc   = bins(1:end-1)+0.5*h;    % <- the bin centers ... 

  vk=histc(x,bins); vk(end)=[];
  fhat = vk/(n*h); % normalize: 

  fh=figure; hold on; 
  ah1=stairs( bins, [fhat;fhat(end)], '-b' ); grid on; 
  xlabel('spatial variable'); ylabel('probability distribution'); 

  % 2) using the normal reference rule:
  % 
  st  = std(x); 
  hNR = (24*(st^(3))*sqrt(pi) / n )^(1/3); 
  h = hNR; 

  nbin = ceil(rng/h); 
  bins = t0:h:(nbin*h+t0);       % <- the bin edges ... 
  bc   = bins(1:end-1)+0.5*h;    % <- the bin centers ... 

  vk=histc(x,bins); vk(end)=[];
  fhat = vk/(n*h); % normalize: 

  figure(fh); 
  ah2=stairs( bins, [fhat;fhat(end)], '-r' ); grid on; 
  xlabel('spatial variable'); ylabel('probability distribution'); 


  % 3) using Scott's rule:
  % 
  st  = std(x); 
  hSc = 3.5 * st * n^(-1/3); 
  h   = hSc; 

  nbin = ceil(rng/h); 
  bins = t0:h:(nbin*h+t0);       % <- the bin edges ... 
  bc   = bins(1:end-1)+0.5*h;    % <- the bin centers ... 

  vk=histc(x,bins); vk(end)=[];
  fhat = vk/(n*h); % normalize: 

  figure(fh); 
  ah3=stairs( bins, [fhat;fhat(end)], '-g' ); grid on; 
  xlabel('spatial variable'); ylabel('probability distribution'); 


  % 4) using the freedman-diaconis rule:
  %

  hFD = 2*iqr(x)*n^(-1/3); 
  h = hFD; 

  nbin = ceil(rng/h); 
  bins = t0:h:(nbin*h+t0);       % <- the bin edges ... 
  bc   = bins(1:end-1)+0.5*h;    % <- the bin centers ... 

  vk=histc(x,bins); vk(end)=[];
  fhat = vk/(n*h); % normalize: 

  figure(fh); 
  ah4=stairs( bins, [fhat;fhat(end)], '-k' ); grid on; 
  xlabel('spatial variable'); ylabel('probability distribution'); 


  legend( [ah1,ah2,ah3,ah4], {['Sturges'' =',num2str(hS)],['Normal Reference ; h=',num2str(hNR)],['Scott''s ; h=',num2str(hS)],['Freedman-Diaconis ; h=',num2str(hFD)]}, 'location', 'northeast' ); 

  if( pi == 1 )
    title('various histogram approximations: galaxy EastWest data');
    axis( [-40 +40 0, 0.055 ] ); 
    saveas(gcf,'../../WriteUp/Graphics/Chapter9/prob_9_3_galaxy','epsc');  
  else
    title('various histogram approximations: forearm data');
    axis( [ 15, 23, 0, 0.5 ] )
    saveas(gcf,'../../WriteUp/Graphics/Chapter9/prob_9_3_forearm','epsc');  
  end
  
end