%-----------------------------------------------------------------------
% Plots input data as a histogram assumes the distribution is Gaussian.
% It also plots a Gaussian distribution for the computed mean and
% standard deviation of the data. Returns the statistics of the data:
%
% Chi  -  Chi Square statistic
% dof  -  Degrees of freedom of the Chi-Square Test
% Prob -  Confidence limit that the data is Gaussian [%]
% mu   -  Mean
% sig  -  Standard Deviation
% Num  -  Number of data points
%
% [Chi, dof, Prob, mu, sig, Num] = % plotHist(data,binSize,minBin, ...
%                                             maxBin,TextOn);
%
% Inputs:
%
% data    - data to be analyzed/plotted
% binSize - width of each bin in the histogram (units are the same as
%           the input data.
% minBin  - Minimum histogram bin, if no argument will find the minimum
%           of the data
% maxBin  - Maximum histogram bin, if no argument and there is a minBin
%           that is negative will take the value abs(minBin). Otherwise
%           if there is no value passed the algorithm will find the
%           maximum of the data.
% TextOn  - Annotates the plot with the number of elemends in each bin
% 
% Stuart Riley
% 7/29/01
%-----------------------------------------------------------------------

function [Chi, dof, Prob, mu, sig, Num] = ...
                            plotHist(data,binSize,minBin,maxBin,TextOn);

TRUE  = 1;
FALSE = 0;

% global ChiTestIn;
% global ChiNorm;
if(nargin==0)
  print 'error no data defined';
  return;
elseif(nargin==1)
  binSize = 1;
  % No range set to the limits of the data.
  minBin  = floor(min(data));
  maxBin  = ceil (max(data));
  TextOn  = TRUE;
elseif(nargin==2)
  % No range set to the limits of the data.
  minBin  = floor(min(data));
  maxBin  = ceil (max(data));
  TextOn  = TRUE;
elseif(nargin==3)
  % Only min defined, set max to be abs(min)
  if(minBin < 0)
    maxBin = -minBin;
  else
    maxBin = ceil( max(data));
  end
  TextOn  = TRUE;
elseif(nargin==3)
  TextOn  = TRUE;
end

nbins = minBin:binSize:maxBin;

[n,xbin]=hist(data,nbins);

Num     = sum(n); % Number of data points
ChiTestIn = n;
% Convert to percentage
n = 100 * n ./ Num;

% Mean of the assumed Normal Distribution
mu  = mean(data); 
% Standard Deviation of the assumed Normal Distribution
sig = std(data); 

H = bar(xbin,n,1); % Plots the histogram.  No gap between bars.
PlotHandlesMat = H;
legendCell = { 'Data Distribution' } ;
hold on

%
% Compute the integral across the bin.
% 
y =   LocalNormCDF(nbins+binSize/2,mu,sig) ...
    - LocalNormCDF(nbins-binSize/2,mu,sig);

% Scale to the number of data points for the Chi-square test
ChiNorm = y .* Num;
% ChiNorm   = y;
% ChiTestIn = ChiTestIn ./ Num;

% Finally conver to percent
y = y * 100; 
% plot(nbins,y,'r','LineWidth',2); % Plots density line over histogram.
[AX,H1,H2] = plotyy(nbins,y,nbins,cumsum(n)); % Plots density line over histogram.

set(get(AX(2),'Ylabel'),'String','Cumulative Distribution [%]')
A = axis(AX(2));
A(1) = minBin;
A(2) = maxBin;
A(4) = 100;
axis(AX(2),A);

A = axis(AX(1));
A(1) = minBin;
A(2) = maxBin;
A(4) = ceil(max(n)) + 1;
axis(AX(1),A);

legendCell = [ legendCell { ['Normal Distribution: \mu=',num2str(mu,'%.1f'),' \sigma=',num2str(sig,'%.1f')]}];
legendCell = [ legendCell { 'Data CDF' } ] ;  % Append legend text.

PlotHandlesMat = [ PlotHandlesMat H1 ] ;
PlotHandlesMat = [ PlotHandlesMat H2 ] ;
set(H1,'color','r');
set(H1,'LineWidth',2);
set(H2,'color','g');
set(H2,'LineWidth',2);


grid
ylabel('Distribution [%]');

legend( PlotHandlesMat( 1, :) , legendCell , 2);

if(TextOn == TRUE)
  for i=1:length(nbins);
    TextArray(i) = {num2str(ChiTestIn(i))};
  end
  text(nbins,n+.05,TextArray,'Rotation',90, ...
                             'FontSize',8, ...
                             'FontWeight', 'bold' ) ;
end

% Compute the Chi-square statistic

TOL = 0.000001;
i = find( ( abs(ChiTestIn) > TOL) | (abs(ChiNorm) > TOL) );
% Chi = sum(  ( ( ChiTestIn(i) - ChiNorm(i) ).^2 )./ ChiNorm(i));
Chi = NaN;
dof = length(i) - 1;

% Chi = sum(  ( ( ChiTestIn - ChiNorm ).^2 )./ ChiNorm);
% dof = length(ChiTestIn) - 1;

% Large Chi indicates that the data is unlikely to be from the
% distribution.

% Compute the likelihood from the Chi^2 statistic
% A High probability means the distributions are likely to be
% from the same distribution.

Prob = NaN; % Chi2cdf is only available in the stats library.
% Prob = 1 - chi2cdf(Chi,dof);


% [H, P, KSTAT, CV] = kstest(data, [nbins; normcdf(nbins,mu,sig)]')
% 1 - gampdf(Chi/2, dof/2);


%
% ------------ Local Functions ---------
%

function dist = LocalNormCDF(x, mu, sig);
  % Normal Cumulative Distribution Function
  dist = (1 + (erf((x-mu)/(sig*sqrt(2)))))./2;

