function [f, chi2] = mynormfit(y, x, c)
% mynormfit.m  1/23/06
% [f, chi2] = mynormfit(y,x,c)
% plots a histogram of the data in the vector y,
% using the vector of class limits given in x.
% The class intervals must all be equal, and the
% range of x should include most of the data.  n is the
% vector of frequencies in each class.  Also plots
% the fitted normal distribution, and calculates ChiSquared
% If c = 1, then the histogram will be color-filled
% written by Gerry Middleton, September 1996.
%tic
ni = length(x) - 1;    % number of classes
N = length(y);         % number of data
dx = (x(2) - x(1))/2;  % find half the class interval
xmp = x + dx;          % determine class mid-points
xmp = xmp(1:ni);       % discard last mid-point
[f,xmp] = hist(y,xmp);
nn = [f 0];
if c == 1              % for filled histogram
   cc = [f;f];         % ff and xx give coordinates of polygon
   cc = reshape(cc,1,2*ni);
   cc = [0 cc 0];
   xx = [x;x];
   xx = reshape(xx,1,(2*ni + 2));
   fill(xx,cc,'g');    % change 'g' to another color if preferred
   line(xx,cc);
   hold on; 
else                   % for unfilled histogram
   stairs(x, nn);
   hold on;
   plot([x(1) x(1)], [0 nn(1)]);
end
hold on;
m = mean(y);           % mean value
s = std(y);            % standard deviation
p = 0.5 + 0.5*erf((x-m)/(sqrt(2)*s));  % cum norm prob of class limits
pt = p(ni+1) - p(1);   % total probability in classes
p2 = p(2:ni+1);
pc = p2 - p(1:ni);     % prob of each class
fc = pc*N*pt;          % theoretical frequency in each class
chi2 = sum((f - fc).^2./fc); % calculate Chi Square
dx2 = (max(x) - min(x))/100; % now plot Normal curve
x2 = [min(x):dx2:max(x)];
x2m = x2(1:100) + dx2/2; % 100 midpoints
p3 = 0.5 + 0.5*erf((x2-m)/(sqrt(2)*s));
p3 = p3(2:101) - p3(1:100);  
fc2 = p3*pt*N*100/ni;
plot(x2m,fc2,'r','linewidth',1);
title('Histogram and fitted Normal curve');
xlabel('x'), ylabel('frequency');         
hold off
%t = toc