function H=chi2test(x, alpha)
%   CHI2TEST: Single sample Pearson Chi Square goodness-of-fit hypothesis test.
%   H=CHI2TEST(X,ALPHA) performs the particular case of Pearson Chi Square
%   test to determine whether the null hypothesis of composite normality PDF is 
%   a reasonable assumption regarding the population distribution of a random sample X
%   with the desired significance level ALPHA.
%
%   H indicates the result of the hypothesis test according to the MATLAB rules 
%   of conditional statements:
%   H=1 => Do not reject the null hypothesis at significance level ALPHA.
%   H=0 => Reject the null hypothesis at significance level ALPHA.
% 
%   The Chi Square hypotheses and test statistic in this particular case are:
% 
%   Null Hypothesis:        X is normal with unknown mean and variance.
%   Alternative Hypothesis: X is not normal.
%
%   The random sample X is shifted by its estimated mean and normalized by its
%   estimated standard deviation. The tested bins XP of the assumed normal distribution
%   are chosen [-inf, -1.6:0.4:1.6, inf] to avoid unsufficient statistics. 
%   
%   Let E(x) be the expected frequency X falls within XP according to the normal
%   distribution and O(x) be the observed frequency. The Pearson statistic, 
%   X2=SUM((E(x)-O(x))^2/E(x)) distributes Chi Square with length(XP)-3 degrees
%   of freedom. 
%
%   The decision to reject the null hypothesis is taken when the P value (probability that Chi2
%   random value with length(XP)-3 degrees of freedomd is greater than X2) is less than 
%   significance level ALPHA. 
%
%   X must be a row vector representing a random sample. ALPHA must be a scalar.
%   The function doesn't check the formats of X and ALPHA, as well as a number of the
%   input and output parameters.
%
%   The asymptotic limit of the Chi Square test presented is reached when LENGTH(X)>90.
%
% Acknowledge: Dr. S. Loyka
%
% Author: G. Levin, May, 2003.
%
% References:
%   W. T. Eadie, D. Drijard, F. E. James, M Roos and B. Sadoulet, "Statistical Methods
%   in Experimental Physics", North-Holland, Sec. Reprint, 1982.

%Normalize x
N=length(x);
x=(x-mean(x))/std(x); %standardization

xp=[-inf, -1.6:.4:1.6, inf]; %tested bins
E=0.5*N*diff(erfc(-xp/sqrt(2))); %expected frequency
S=histc(x, xp); 
O=S(1:end-1); %%observed frequency
%plot(xp(2:end),E,'k-',xp(2:end),O,'k.');
x2=sum((E-O).^2./E); %statistics

pval=1-gammainc(x2/2,(length(O)-3)/2); %p value

H=(pval>=alpha);