function [table, chi2, p, labels] = crosstab(varargin)
%CROSSTAB Cross-tabulation of column vectors.
%   CROSSTAB(COL1,COL2,...) takes vectors, string arrays, or cell
%   arrays of strings, and returns an array, TABLE, of crosstabs.
%   Element (i,j,...) of TABLE contains the count of all instances
%   where COL1 = i, COL2 = j, and so on.
%
%   [TABLE, CHI2, P] = CROSSTAB(...) also returns the chisquare
%   statistic, CHI2, for testing independence of each dimension of
%   TABLE. (That is, that the proportion of items falling in any
%   cell is equal to the product of the proportion in that row,
%   times the proportion in that column, and so on.)  The scalar,
%   P, is the significance level of the test.  Values of P near
%   zero cast doubt on the assumption of independence.
%
%   [TABLE, CHI2, P, LABELS] = CROSSTAB(...) also returns a cell
%   array of labels for TABLE.  The entries in the first column of
%   LABELS are labels for the rows of TABLE, the entries in the
%   second colum are labels for the columns, and so on.

%   B.A. Jones 3-5-95
%   Extensive re-write 9-23-99 by Tom Lane
%   Copyright 1993-2000 The MathWorks, Inc. 
%   $Revision: 2.9 $  $Date: 2000/05/26 17:28:41 $

sz = zeros(1,nargin);
if (nargout > 3), labels = cell(0, nargin); end
nonan = [];
M = [];
for j=1:nargin
   [g1,g2] = grp2idx(varargin{j});
   ng = size(g2,1);
   if (nargout > 3), labels(1:ng,j) = g2; end
   sz(j) = ng;
   n = length(g1);
   if (j == 1)
      n1 = n;
      nonan = ~isnan(g1);
      M = zeros(n,nargin);
      M(:,1) = g1;
   elseif (n ~= n1)
      error(sprintf('Arguments 1 and %d have different lengths.', j));
   else
      nonan = nonan & ~isnan(g1);
      M(:,j) = g1;
   end
end

M = M(nonan,:);
n = size(M,1);

if (length(sz) > 1)
   table = zeros(sz);
else
   table = zeros(sz, 1);
end

for k = 1:n
   N = num2cell(M(k,:));
   table(N{:}) = table(N{:}) + 1;
end
 
if (nargout > 1 & nargin > 1)
   expected = zeros(size(table));
   expected(:) = n;
   T = table;
   szv = sz;
   permv = [(2:nargin) 1];
   for j = 1:nargin
      sz1 = szv(1);
      sz2 = prod(szv(2:end));
      frac = sum(reshape(T, sz1, sz2), 2) / n;
      frac = reshape(repmat(frac,1,sz2),szv);
      expected = expected .* frac;
      expected = shiftdim(expected, 1);
      szv = szv(permv);
      T = shiftdim(T, 1);
   end

   chi2 = (table - expected).^ 2 ./ expected;
   chi2 = sum(chi2(:));
   df = prod(sz) - (1 + sum(sz - 1));
   p = 1 - chi2cdf(chi2,df);

elseif (nargout > 1)      % nargin <= 1
   chi2 = NaN;
   p = NaN;
end