%% 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 Z (scostamenti standardizzati)
% Z= 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 Z 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
% Cancellare tutte le variabili in memoria.
% Rileggere di nuovo i dati dal file di Excel e calcolare l'analisi delle
% corrispondenze utilizzando la funzione CorAna del toolbox FSDA.
%

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

X = readtable('corrisp.xlsx', 'Sheet','dati','range','A1:D427','ReadRowNames',true);

%% Creo la tabella di contingenza tra le variabili professione e tipo di acquisto

[N,chi2,pvalchi2,labels] =crosstab(X.Professione,X.Tipo_di_Acquisto);
% N è la la matrice che contiene le frequenze per ogni combinazione di
% Professione e Tipo_di_Acquisto (tabella di contingenza)
% chi2 è il valore del test chi quadrato
% pvalchi2 è il valore del p-value del test chi quadrato di
% label è un array di dimensione 13 x 4 che contiene nella prima colonne le
% etichette delle modalità di riga (Professione) e nella seconda colonna
% (nelle prime 4 posizioni) le etichette delle modalità di colonna (Tipo
% di acquisto)

[I,J]=size(N);

NamesRows = labels(1:I,1);
NamesCols = labels(1:J,2);

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


%% 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";
Ntable=Ntable(:, {'Last minute', 'Giornaliero', 'Ricorrenze', 'Settimanale'});


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


%% 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 = N./sum(N,2);

% ProfilesCols = matrice che contiene i profili colonna
% Si divide ogni colonna per il suo totale
ProfilesCols = N./sum(N,1);


%% 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)

% r= vettore che contiene le masse di righe = centroidi dei profili colonna
r=sum(N,2)/n;

% Masse di colonna = centroidi dei profili di riga
c=(sum(N,1)/n)';

%% 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);


%% 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
cchk=(ProfilesRows')*r;
rchk=ProfilesCols*c;

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

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

%% Effettuare la scomposizione in valori singolari della matrice Z 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 dif Z
[U,Gam,V] = svd(Z,'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);


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

% TotalInertia = sum_i sum_j (pij - ricj)^2 / ricj = X^2/n
TotalInertia     = sum(sum(Z.^2));

Chi2stat     = n*TotalInertia;

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

% 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 = varianza totale della tabella
%   di contingenza)
% - 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);


%% 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;

%% 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 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={'Grafico di analisi delle corrispondenze.' , ...
    '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')
xline(0)
yline(0)

%% 
% Cancellare tutte le variabili in memoria.
% Rileggere di nuovo i dati dal file di Excel e calcolare l'analisi delle
% corrispondenze utilizzando la funzione CorAna del toolbox FSDA.

clear
X = readtable('corrisp.xlsx', 'Sheet','dati','range','A1:D427','ReadRowNames',true);
out=CorAna(X,'datamatrix',true);