function [stats] = bowker( data )
% Test for Internal Symmetry
%Bowker, A. H. (1948).
%A test for symmetry in contingency tables. 
%J. Am. Stat. Assoc, 43, 572 - 574.
%  
% Example of use
%    eyes=[...
%            1520   266     124      66  ; ... 
%            234     1512    432      78  ;... 
%            117     362  1772     205 ; ...
%             36      82    179      492];
%  [stats] = bowker(eyes)
%     
%     
[k,l] = size(data);
if (k ~= l) 
     error('Data Matrix Needs to be Square') 
end
n=sum(sum(data));
bowk=0; mjs=0; bowkcorrected=0;
k1 = k-1;
for i = 1:k1
    i1=i+1;
    for j = i1: k 
        %Bowker Stat
        bowij = (data(i,j)-data(j,i))^2 /(data(i,j)+data(j,i)+eps);            
        bowk = bowk + bowij;
        % Bowker Corrected (Edwards 1948 Correction)
        bowijc = (abs(data(i,j)-data(j,i))-(abs(data(i,j)+data(j,i))>0))^2 /...
                (data(i,j)+data(j,i)+eps);            
        bowkcorrected = bowkcorrected + bowijc;
        % May-Johnson Stat 
        mj=n*(data(i,j)-data(j,i))^2 / ...
        (n*(data(i,j)+data(j,i)+eps)-(data(i,j)-data(j,i))^2 );  
        mjs = mjs + mj;
    end
end
df = 1/2 * k * (k-1); %common df
% Bowker
pvalb = 1 - chi2cdf(bowk, df);
stats.bowkerchi2 = bowk;
stats.bowkerpval = pvalb;
% Bowker Corrected
pvalbcorrected = 1 - chi2cdf(bowkcorrected, df);
stats.bowkercorrectedchi2 = bowkcorrected;
stats.bowkercorrectedpval = pvalbcorrected;
% May-Johnson
pvalmj = 1 - chi2cdf(mjs, df);
stats.mayjohnsonchi2 = mjs;
stats.mayjohnsonpval = pvalmj;
end
% Bowker, A.H. (1948): A Test for symmetry in Contingency Tables. Journal
% American Statistical Association, 43, 572-574.

% May, W.L., Johnson, W.D. (2001): Symmetry in Square Contingency Tables:
% Tests of Hypotheses and Confidence Interval Construction. 
% Journal of Biopharmaceutical Statistics, 11, 23-33.

% Edwards, A.  L. (1948). Note on the "Correction for Continuity "in Testing the
% Significance of the Difference Between Correlated Proportions. Psychometrika, 13,
% 185-187.

