function [chi2, df, pvalue, expected] = tablerxcxp(observed)
% % 3D Contingency Table r x c x p for testing the independence
% %
% % Input:
% %  observed - r x c x p 3D array of observed frequencies
% %  
% % Output
% %    chi2 - statistics (approx distributed as chi-square with
% %                       (r-1)(c-1) degrees of freedom)
% %    df - degrees of freedom for chi^2 under H_0
% %    pvalue  - p - value
% %    expected -  3D array  of expected frequencies
% % Example of use: 
% obs= cat(3, [44 12; 28 22],[38 10; 20 18]);
% [chi2,df, pvalue,expected]=tablerxcxp(obs)
% %------------------------------------------------------------

[r c p] = size(observed);
    rowsums = sum(sum(permute(observed, [3 2 1])));
    colsums = sum(sum(permute(observed, [3 1 2])));
    pagsums = sum(sum(observed));
    n = sum(pagsums);
 for i = 1:r
     for j = 1:c
         for k=1:p
             expected(i,j,k)= rowsums(i)*colsums(j)*pagsums(k)/n^2;
         end
     end
 end
 chi2 = sum(sum(sum( (observed - expected).^2 ./expected )));
 df = r*c*p - r -c - p + 2; 
 pvalue = 1-chi2cdf(chi2, df);
%-------------------------------------------------------------------
% Copyright (C) 2009 by BVidakovic
% Ver 0.12   8/7/2009% 

