% mcatrad: Multiple Correspondence analysis of an Indicator Matrix
% 
% Input:
% D - Original data matrix of the form "objects x variables"

% The standard MCA algorithm is applied to the corresponding indicator matrix of
% the form "objects x variable categories" 
% The algorithm results in the standard and principal coordinates  
% of the row and column profiles.

% [Sr,Sc,Pr,Pc] = mcastd(D);

% Output:
% Sr - Row standard coordinates
% Sc - Column standard coordinates
% Pr - Row principal coordinates 
% Pc - Column principal coordinates

% -- Example --

% D = [2     2     1     
%      3     2     3     
%      2     1     2     
%      3     1     2     
%      3     1     2     
%      2     2     3   
%      2     2     2    
%      2     1     1     
%      3     2     1   
%      2     2     1   
%      3     2     3   
%      1     1     3   
%      1     1     3   
%      1     2     3   
%      1     2     2   
%      2     2     2]

% Author: Angelos Markos 

function [Sr,Sc,Pr,Pc] = mcastd(D)

%Calculate Indicator Matrix
[m,p] = size(D);
% variable levels
nlev = [];
for index = 1:p
    nlev = [nlev,length(unique(D(:,index)))];
end
% total number of categories
q = sum(nlev);
% cumsum of cats
offset = [0,cumsum(nlev)];
Z = zeros(m,q);
for i = 1:m
    for j = 1:p
        Z(i,(offset(j) + D(i,j))) = 1;
    end
end

[r,c] = size(Z);

% Grand total of counts and marginal totals of probabilities
total = sum(sum(Z));                
coltot = sum(Z);                    
rowtot = sum(Z,2)';

% Step 1: Calculate the matrix Rz of standardized residuals
Rz = (Z./sqrt(rowtot'*coltot))-(sqrt(rowtot'*coltot)/total);
clear Z;

% Step 2: Calculate the SVD of Rz (thin form)  
if r>=c 
    [a,w,v] = svd(Rz,0); 
else
    [v,w,a] = svd(Rz',0);
end

% Maximum number of non trivial axes 
np = min(c - p, r-1);
w=w(1:np,1:np);

% Step 3: Standard coordinates of rows and columns
 %row and column masses
ri = rowtot / total;
cj = coltot / total;     
Sr = a(:,1:np)./(sqrt(ri)'*ones(1,np));
Sc = v(:,1:np)./(sqrt(cj)'*ones(1,np));

% Step 4: Principal coordinates of rows and columns
Pr = Sr*w;
Pc = Sc*w;