function [stats] = stuartmaxwell( data )
% Test for Marginal Symmetry
% Stuart-Maxwell test and assymptotically
% equivalent Bhapkar test.
% Example of use
%    eyes=[...
%            1520   266     124      66  ; ... 
%            234     1512    432      78  ;... 
%            117     362  1772     205 ; ...
%             36      82    179      492];
%  [stats] = stuartmaxwell (eyes)
%     
%     
[k,l] = size(data);
if (k ~= l) 
     error('Data Matrix Needs to be Square') 
end
k1 = k-1;
ndoti = sum(data);
nidot = sum(data');
n=sum(sum(data));
di = ndoti - nidot;
di1 = di(1:k1);
for i = 1:k1
    for j = 1: k1
        if(i == j)
            si(i,j) = nidot(i) + ndoti(i) - 2 * data(i,i);            
        else
            si(i,j) = -(data(i,j) + data(j,i));
        end
    end
end
chi2 = di1 * inv(si) * di1';
bsi= si - 1/n * di1' * di1;
bchi2 = di1 * inv(bsi) * di1' ;
  % or bchi2=chi2/(1-chi2/n) directly
pval = 1 - chi2cdf(chi2, k-1);
bpval = 1 - chi2cdf(bchi2, k-1);
stats.stuartmaxwellchi2 = chi2;
stats.stuartmaxwellpval = pval;
stats.bhapkarchi2= bchi2;
stats.bhapkarpval=bpval;
end