#
#   function [C,d,M,mm]=burt(H)
#
#   Computes the bivariate marginals matrix (Burt's matrix)
#   for a data matrix.
#
#   Input: H, data matrix, with n rows and m columns.
#          It is assumed that each variable (= column of H)
#          is coded in integer values, with minimum = 1
#          and maximum = number of category levels.
#
#   Output:  C, the bivariate marginals matrix (Burt's matrix)
#               C = G'*G, where G is the complete indicator
#               matrix. G is implicitely computed but not saved.
#            d, the univariate marginals vector (row vector
#               whose entries are the diagonal entries in C
#               = sum of columns of the indicator matrix G).
#            M, row vector (1,m). Each entry of M equals the
#               number of category levels of the corresponding
#               column of H.
#            mm = sum(M), total number of category levels
#               (C and D have mm rows and mm columns).
#
function [C,d,M,mm]=burt(H)
   [n,m]=size(H);
   M=max(H);
   mm=sum(M);
   Gi=zeros(1,mm);
   C=zeros(mm,mm);
   d=zeros(1,mm);
   for i=1:n,
      jmax=0;
      for alpha=1:m,
         jmin=jmax+1;
         jmax=jmin+M(alpha)-1;
         tt=zeros(1,M(alpha));
         tt(H(i,alpha))=1;
         Gi(jmin:jmax)=tt;
      end
      C=C+Gi'*Gi;
      d=d+Gi;
   end
endfunction