%-----------------------------------------------------------------------------
%
% $Id: chisq.m,v 1.1 1999/09/22 07:34:21 philipc Exp $
%
% Program to perform chi^2 test 
%
%-----------------------------------------------------------------------------
function chi=chisq(x, n, mu, sigma)

numBins = size(x, 2);

% Add up how many points there are in the histogram

N = sum(n);

%
% For each frequency, we need the probability that a r.v. falls into a
% given bin. We obtain this by integrating the Gaussian d.f. over the
% range of the bin. The centers of the bins are given by x(1) to x(numBins),
% so the ranges for the bins are given by
%
% (-Inf, (x(1)+x(2))/2] ((x(1)+x(2))/2, (x(2)+x(3))/2] . . .
%        . . . ((x(k-1) + x(k))/2, +Inf)
%
% noting that the two end bins need special treatment because they extend to
% infinity.
%
% The integral is found by integrating the Gaussian d.f. with zero mean and
% unit s.d., and then scaling appropriately, using the Matlab error function
% erf().
%

% First calculate the upper and lower ends of the bins, and translate them
% to zero mean and unit s.d.

a = zeros(1, numBins-1);
for i = 1:numBins-1
    a(i) = (0.5*(x(i) + x(i+1)) - mu)/(sigma*sqrt(2));
end;

% Now calculate the areas (probabilities) for interior bins

prob = zeros(1, numBins);

for i = 2:numBins-1
    prob(i) = 0.5*(erf(a(i)) - erf(a(i-1)));
end;

% Now the last bin from a(end) to +Inf

prob(end) = 0.5*erfc(a(end));

% Now the first bin from a(1) to -Inf (actually, from -a(1) to +Inf)

prob(1) = 0.5*erfc(-a(1));

% Now calculate the chi^2 statistic

chicomp = zeros(1, numBins);

for i = 1:numBins
   chicomp(i) = (n(i) - N*prob(i))^2/(N*prob(i));
end;

chi = sum(chicomp);

%
% end function chisq
%