%% TESTO
% Importare dentro MATLAB il dataset contenuto nel foglio dati del file
% corrisp.xlsx
% Creare la tabella di contingenza tra le variabili
% Professione (sulle righe) e Tipo_di_Acquisto (sulle colonne).
% Riodinare le modalità di riga della tabella di contingenza in modo da
% inserire la modalità "Altro" come ultima riga. Riordinare anche le
% relative labels.
% Riodinare le modalità di colonna della tabella di contingenza in modo da
% seguire il seguente ordine:
% "last minute", "giornaliero", "ricorrenze", "settimanale;
% Per maggiore leggibilità trasformare la matrice contenente la tabella di
% contingenza in formato table.
% Calcolare la matrice delle corrispondenze P che contiene le frequenze
% relative.
% Calcolare i vettori r e c che contengono rispettivamente le masse di riga
% e di colonna (tramite somme oppure tramite moltiplicazioni matriciali)
% Costruire le due matrici Dr e Dc che contengono rispettivamente i profili
% medi di colonna e riga sulla diagonale principale.
% Calcolare le matrici dei profili riga e dei profili colonna (tramite
% moltiplicazioni matriciali).
% Verificare che le masse di colonna non altro che i profili medi di riga e
% che le masse di riga non sono altro che i profili medi di colonna.
% Costruire la matrice S (scostamenti standardizzati)
% S= Dr^(-1/2) * (P - r * c') * Dc^(-1/2)
% Partendo dalla matrice S calcolare l'inerzia totale discutendo la sua
% relazione con l'indice di Cramer
% Mostrare che l'inerzia totale si può ottenere anche come somma
% ponderata delle distanze di ciascun profilo riga dal profilo medio
% oppure come somma ponderata di ciascun profilo colonna dal profilo medio.
% Effettuare la scomposizione in valori singolari della matrice S e trovare
% le coordinate dei punti riga e colonna da rappresentare nel grafico
% Dare un giudizio sulla bontà dell'analisi calcolando il contributo della
% generica componente i-esima all'inerzia totale.
% Calcolare gli score di riga e di colonna delle prime 3 dimensioni latenti.
% Calcolare i contributi di ogni punto (riga o colonna) all'inerzia della
% dimensione e i contributi della dimensioni all'inerzia dei punti riga e
% colonna
% Rappresentare in un diagramma a dispersione gli score di riga e quelli di
% colonna utilizzando simboli diversi.
% Discutere i risultati ottenuti

%% -----------------------------------------
% Soluzione
%-----------------------------------------

%% Importare dentro MATLAB il dataset contenuto nel foglio dati del file
% corrisp.xlsx
[~, X] = xlsread('corrisp.xlsx','dati','B2:D427');


%% Creare la tabella di contingenza tra le variabili
% Professione (sulle righe) e Tipo_di_Acquisto (sulle colonne).
[N,~,~,labels] =crosstab(X(:,1),X(:,2));
[I,J]=size(N);

%% Riodinare le modalità di riga della tabella di contingenza in modo da
% inserire la modalità "Altro" come ultima riga. Riordinare anche le
% relative labels.
% Riodinare le modalità di colonna della tabella di contingenza in modo da
% seguire il seguente ordine:
% "last minute", "giornaliero", "ricorrenze", "settimanale;

% La modalità altro altra viene spostata come ultima riga
N=[N(2:end,:); N(1,:)];
labels(:,1)=[labels(2:end,1); labels(1,1)];

% Le colonne vengono inserite con ordine
% last minute, giornaliero ricorr. settimanale
neworder=[3 1 4 2];
N=N(:,neworder);
labels(1:4,2)=labels(neworder,2);

%% Per maggiore leggibilità trasformare la matrice contenente la tabella di
% contingenza in formato table.

% Le due istruzioni che seguno non sono più necessarie in MATLAB 2019b
NamesRows = labels(1:I,1);
NamesCols = labels(1:J,2);

% Dall'array N alla table Ntable
Ntable=array2table(N,'RowNames',NamesRows,'VariableNames',NamesCols);
disp(Ntable)


%% Calcolare la matrice delle corrispondenze P che contiene le frequenze
% relative.

% n= numero di unità del campione
n=sum(sum(N));

% P = (matrice di corrispondenza =correspondence matrix) contiene le
% frequenze relative f_ij
P = (1/n) * N;

%% Calcolare i vettori r e c che contengono rispettivamente le masse di riga
% e di colonna (tramite somme oppure tramite moltiplicazioni matriciali)

%vectors of ones
onesI = ones(I,1);
onesJ = ones(J,1);

% r= vector which contains row masses = centroids of the column profiles
r  = P * onesJ ;
% senza la moltiplicazione matriciale
rchk=sum(P,2);


% Masse di colonna = centroidi dei profili di riga r' * Dr^(-1) * P = 1' * P = c'
c  = (onesI' * P)';
% senza la moltiplicazione matriciale
cchk=sum(P,1)';

%% Costruire le due matrici Dr e Dc
% Queste matrici contengono rispettivamente i profili
% medi di colonna e riga sulla diagonale principale.
Dr = diag(r);
Dc = diag(c);



%% Calcolare le matrici dei profili riga e dei profili colonna (tramite
% moltiplicazioni matriciali).

% ProfilesRows = matrice che contiene i profili riga
% Si divide ogni riga per il suo totale
ProfilesRows = Dr^(-1) * P;
% Modo alternativo per trovare la matrice dei profili riga
% P./sum(P,2)
% Modo alternativo per trovare la matrice dei profili riga
% P./r

% ProfilesCols = matrice che contiene i profili colonna
% Si divide ogni colonna per il suo totale
ProfilesCols = P * Dc^(-1);
% Modo alternativo per trovare la matrice dei profili colonna
% P./c'

%% Verificare che le masse di colonna non altro che i profili medi di riga e
% che le masse di riga non sono altro che i profili medi di colonna
cchk1=ProfilesRows'*r;
cchk2= sum(ProfilesRows.*r,1)';

rchk1=ProfilesCols*c;
rchk2=sum(ProfilesCols.*c',2);

%% Costruire la matrice S (scostamenti standardizzati)
% Nhat = freuenze attese sotto l'ipotesi di indipendenza
Nhat=(r*c')*n;

% S   = Dr^(1/2) * (Dr^(-1) * P - r * c') * Dc^(-1/2) = Dr^(-1/2) * (P - r * c') * Dc^(-1/2);
% sij = sqrt( (p_{ij} - r_ic_j)^2 /r_ic_j ) =(p_{ij} - r_ic_j)/sqrt(r_ic_j)
% Chi-square distances
S     =  Dr^(-1/2) * (P - r * c') * Dc^(-1/2);


%% Partendo dalla matrice S calcolare l'inerzia totale discutendo la sua
% relazione con l'indice di Cramer

% Total inertia = somma dei quadrati degli elementi della matrice S =
% varianza totale della tabella di contingenza
% TotalInertia = sum_i sum_j (pij - ricj)^2 / ricj = X^2/n
TotalInertia     = sum(sum(S.^2));

Chi2stat     = n*TotalInertia;

% Cramer's V
CramerV = sqrt(Chi2stat/(n*(min(I,J)-1)));

%% Inerzia ottenuta come somma ponderata delle distanze di ciascun profilo
% dal profilo medio

% distI = vettore che contiene nell'elemento i-esimo la distanza della
% modalità i dal profilo medio
distI=zeros(I,1);

for i=1:I
    distI(i)=sqrt(sum(  (ProfilesRows(i,:)-c').^2./(c')  ));
end
TotalInertiachk=(distI.^2)'*r;


%% Effettuare la scomposizione in valori singolari della matrice S e trovare
% le coordinate dei punti riga e colonna da rappresentare nel grafico
% Dare un giudizio sulla bontà dell'analisi calcolando il contributo della
% generica componente i-esima all'inerzia totale

% SVD of S
[U,Gam,V] = svd(S,'econ');

% K = maximun number of dimensions
K = min(I-1,J-1);

Gam = Gam(1:K,1:K);
U   = U(:,1:K);
V   = V(:,1:K);


% cumsumTotalInertia = cumulative proportion of explained inertia
Gam2 = Gam.^2;
cumsumTotalInertia = cumsum(diag(Gam2))/TotalInertia;

% InertiaExplained è una matrice con quattro colonne.
% - La prima colonna contiene i valori singolari (la somma dei quadrati
%   dei valori singolari è l'inerzia totale)
% - La seconda colonna contiene gli autovalori (ossia i quadrati dei valori singolari).
% (la somma degli autovalori è l'inerzia totale)
% - La terza colonna contiene la varianza spiegata da ciascuna dimensione latente
% - La quarta colonna contiene la varianza cumulata spiegata da ciascuna dimensione
InertiaExplained=[diag(Gam) diag(Gam2) diag(Gam2 / TotalInertia) cumsumTotalInertia];
ColNamesSummary={'Valori_singolari' 'Autovalori' 'Var_spiegata' 'Cum_Var_spiegata'};
RowNamesSummary=strcat(cellstr(repmat('dim_',K,1)), cellstr(num2str((1:K)')));
RowNamesSummary=regexprep(RowNamesSummary,' ','');

InertiaExplainedtable=array2table(InertiaExplained,'VariableNames',ColNamesSummary, ...);
    'RowNames',RowNamesSummary);



%% Calcolare gli score di riga e di colonna delle prime 3 dimensioni latenti.
% e i contributi di ogni punto (riga o colonna) all'inerzia della

% Principal coordinates of rows (scores di riga)
RowsPri     = Dr^(-1/2) * U*Gam;
% Osservazione RowsPri'*Dr*RowsPri = matrice degli autovalori = quadrati
% dei valori singolari sulla diagonale principale
% Contributi delle diverse righe alla determinazione del primo autovalore
disp(RowsPri(:,1).^2.*r/(Gam(1,1)^2))


% Osservazione: la quantità
% sum(RowsPri(:,1).^2.*r)
% è la varianza ponderata della prima dimensione (primo autovalore =Gam(1,1)^2)


% Contributi delle diverse righe alla determinazione del secondo autovalore
disp(RowsPri(:,2).^2.*r/(Gam(2,2)^2))
% Osservazione: la quantità
% sum(RowsPri(:,2).^2.*r)
% è la varianza ponderata della seconda dimensione (secondo autovalore =Gam(2,2)^2)

% Principal coordinates of columns G= (Dc^(-1/2)*V*Gamma)
% Scores di colonna
ColsPri     = Dc^(-1/2) * V*Gam;
% Osservazione ColsPri'*Dc*ColsPri = matrice degli autovalori = quadrati
% dei valori singolari sulla diagonale principale

% Contributi delle diverse colonne alla determinazione del primo autovalore
disp(ColsPri(:,1).^2.*c/(Gam(1,1)^2))
% Osservazione: la quantità
% sum(ColsPri(:,1).^2.*c)
% è la varianza ponderata della prima dimensione (primo autovalore =Gam(1,1)^2)

% Contributi delle diverse colonne alla determinazione del secondo autovalore
disp(ColsPri(:,2).^2.*c/(Gam(2,2)^2))
% Osservazione: la quantità
% sum(ColsPri(:,2).^2.*c)
% è la varianza ponderata della seconda dimensione (secondo autovalore =Gam(2,2)^2)

%% Calcolare i contributi della dimensioni all'inerzia dei punti

% Osservazione: il vettore distI.^2  contiene la distanza al quadrato di
% ogni profilo riga dal profilo medio (inerzia ossia varianza di ogni punto)
% Il vettore
% Contributi della dimensione 1 alla spiegazione dell'inerzia di ogni punto
disp(RowsPri(:,1).^2./distI.^2)
% Contributi della dimensione 2 alla spiegazione dell'inerzia di ogni punto
disp(RowsPri(:,2).^2./distI.^2)




%% Rappresentare in un diagramma a dispersione gli score di riga e quelli di
% colonna utilizzando simboli diversi.
% Discutere i risultati ottenuti
close all
titl={'French symmetrical model: rows and cols in principal coordinates.' , ...
    'Plot of $X=D_r^{-1/2}U \Gamma$ and $Y= D_r^{-1/2} V \Gamma$'};
symbolrows='o';
symbolcols='^';
% Color for symbols and text for rows points
colorrows='b';
% Color for symbols and text for  column rows
colorcols='r';
MarkerSize=14;
hold('on')
plot(RowsPri(:,1),RowsPri(:,2),'LineStyle','none','Marker',symbolrows,'Color', colorrows,'MarkerSize',MarkerSize)
plot(ColsPri(:,1),ColsPri(:,2),'LineStyle','none','Marker',symbolcols,'Color', colorcols,'MarkerSize',MarkerSize)
Lr=labels(1:I,1);
Lc=labels(1:J,2);
dx=0.05;
text(RowsPri(:,1),RowsPri(:,2)+dx,Lr)
text(ColsPri(:,1),ColsPri(:,2)+dx,Lc)

FontName='Times';
FontSizeAxisLabels=12;
title(titl,'Interpreter','Latex');

% Labels for axes
xlabel(['Dimension ',sprintf('%2.0f',1),' (',sprintf('%5.1f',InertiaExplained(1,3)*100),'%)'],'FontName', FontName, 'FontSize', FontSizeAxisLabels);
ylabel(['Dimension ',sprintf('%2.0f',2),' (',sprintf('%5.1f',InertiaExplained(2,3)*100),'%)'],'FontName', FontName, 'FontSize', FontSizeAxisLabels);
axis(gca,'equal')
vv=axis;
line([vv(1);vv(2)],[0;0])
line([0;0],[vv(3);vv(4)])