%% Caricamento dati
clear 
close all
Xca=readtable('esercitazione2020.xlsx','Range','A1:B194');

%% Distribuzione di frequenze variabile Posizione e grafico a barre
% Modo 1
% Trasformo la prima colonna di Xca in variabile categorica e applico la
% funzione summary e la funzione histogram
X1cat=categorical(Xca{:,1});
summary(X1cat);
histogram(X1cat)

% Modo 2
% Applico la funzione tabulate (a1 diventa un cellarray)
a1=tabulate(Xca{:,1});
% La seconda colonna di a1 viene trasformata in double
% freq è un vettore di doubles.
freq=cell2mat(a1(:,2));
bar(freq)
% Aggiungo le etichette sull'asse x (e le ruoto)
set(gca,'xtick',1:length(freq),'xticklabel',a1(:,1),'XTickLabelRotation',60)

% Modo 3
% Applico la funzione tabulate
a1=tabulate(Xca{:,1});
% La seconda colonna di a1 viene trasformata in double
freq=cell2mat(a1(:,2));
% L'istruzione categorical(a1(:,1),a1(:,1)) significa
% trasforma in categorico il vettore a(:,1) e prendi come ordine delle
% modalità da rappresentare nel grafico l'ordine specificato da a1(:,1))
xlabels=categorical(a1(:,1),a1(:,1));
bar(xlabels,freq)

%% NON NELL'ESAME
% esploro la differenza tra l'istruzione 
% categorical(a1(:,1))
% e 
% categorical(a1(:,1),a1(:,1))
subplot(2,1,1)
xlabels=categorical(a1(:,1));
bar(xlabels,freq)
title('L''ordine delle barre è associato all''ordine alfabetico di xlabels')
subplot(2,1,2)

xlabels=categorical(a1(:,1),a1(:,1));
bar(xlabels,freq)
title('L''ordine delle barre è associato all''ordine specificato dentro a1(:,1)')



%% Calcolo indici Cramer e tauyx
Indici=corrNominal(Xca(:,1:2),'datamatrix',true);
disp('Indice di Cramer')
disp(Indici.CramerV)
disp(['Indice tauyx è:' num2str(Indici.tauyx)])
disp(Indici.ConfLimtable('tauyx',:))


%% Analisi delle corrispondenze
Acor=CorAna(Xca(:,1:2),'datamatrix',true');
CorAnaplot(Acor)
% interpretazione dim 1 (da non fumatore a fumatore incallito)
% ==> propensione verso il fumo
% interpretazione dim 2 (manager in basso impiegati in alto)
% ==> Relazione inversa con il grado di importanza
% Inerzia spiegata vicino al 100 per cento. Rappresentazione ottima.

% Mostro la distanza al quadrato moltiplicata per la massa per la riga dei
% Senior Managers ed il profilo medio
% Overview del profilo riga dei Senior Managers
Acor.OverviewRows('Senior_Managers',:)
disp('Distanza al quadrato ponderata per la massa per i SM')
Acor.OverviewRows('Senior_Managers','Inertia')

% Implementazione manuale di questa quantità
% Osservazione il profilo medio di riga è dato dalle masse di colonna
ProfiloRigaSM=Acor.ProfilesRows(1,:);
% ProfiloRigaMedio (come vettore riga) = vettore delle masse di colonna
% trasposto
ProfiloRigaMedio=Acor.c';

Distanza2=sum((ProfiloRigaSM-ProfiloRigaMedio).^2./ProfiloRigaMedio);
MassaRigaSM=Acor.OverviewRows{1,1};
% InerziaSM = inerzia dei senior manager = distanza euclidea al quadrato
% dal profilo medio ponderata con le masse di colonna per la rispettiva massa di riga
% Si può interepretare come (x_i-M)^2 *f_i
InerziaSM=Distanza2*MassaRigaSM;
% La somma delle inserzie dei diverse punti riga è l'inerzia totale ossia
% la variabilità totale della tabella di contingenza

tabCont=Acor.Ntable;
tabCont=tabCont([4 5 3 2 1],:);
corrOrdinal(tabCont)
% L'intervallo di confidenza di gamma contiene il valore zero==>
% associazione tra le due variabili ordinali non significativa

%% Calcolo indice gamma in maniera manuale
% Estrai il numero di righe e di colonne delle tabella di input N
N=table2array(tabCont);
[I,J] = size(N);

% Trovo le coppie discordanti e concordanti 
    C=0;
    D=0;
    for i=1:I
        for j=1:J
            if i<I && j<J
                xsel=N(i+1:I,j+1:J);
                C=C+N(i,j)*sum(xsel(:));
            end
            if j>1 && i<I
                xsel=N(i+1:I,1:j-1);
                D=D+N(i,j)*sum(xsel(:));
            end
        end
    end

% Gamma index
GammaM = (C - D)/(C + D);

%% Caricamento dati
Xpca=readtable('esercitazione2020.xlsx','Range','Q2:V22',...
    'ReadRowNames',true,'Format','auto');
[Z,medie,sigma]=zscore(Xpca{:,:});
figure
boxplot(Z,'Labels',Xpca.Properties.VariableNames)
biplotFS(Xpca)
figure
outPCA=pcaFS(Xpca);

%% Matrice dei diagrammmi di dispersione
spmplot(Xpca)
cor=corr(Xpca{:,:});
array2table(cor,'RowNames',Xpca.Properties.VariableNames,...
    'VariableNames',Xpca.Properties.VariableNames)

%% Analisi con 3 gruppi (metodo delle k medie)
[classif,Cent]=kmeans(Z,3,'Replicates',100);
% Osservazione. Cent contiene già i valori dei centroidi nello spazio delle
% variabili standardizzate

% Y = score non standardizzati
Y=outPCA.score;
etirighe=Xpca.Properties.RowNames;

% Proietto i centroidi nello spazio delle prime due componenti principali
% In altri termini applico ai centroidi le stesse trasformazioni che ho
% applicato alla matrice X per ottenere le componenti principali
% ossia standardizzo moltiplico per i primi due autovettori
centroidiPCA=Cent*outPCA.coeff;

%% Rappresentazione degli score con etichette delle province e centroidi dei gruppi
close all
plot(Y(:,1),Y(:,2),'o')
axislim=axis;
text(Y(:,1),Y(:,2),etirighe)
% Aggiungo l'asse x
line([axislim(1);axislim(2)], [0;0], 'Color','black');
% Aggiungo l'asse y
line([0;0],[axislim(3);axislim(4)], 'Color','black');

hold('on')
plot(centroidiPCA(:,1),centroidiPCA(:,2),'s','MarkerFaceColor','r','LineWidth',5)


%% Rappresentazione degli score + centroidi + gruppo di appartenenza
close all
plot(Y(:,1),Y(:,2),'o')
axislim=axis;
xlabel('Prima componente principale (Indice di povertà )');
ylabel('Seconda componente principale (Indice di malessere delle aziende)');
text(Y(:,1),Y(:,2),num2str(classif))
% Aggiungo l'asse x
line([axislim(1);axislim(2)], [0;0], 'Color','black');
% Aggiungo l'asse y
line([0;0],[axislim(3);axislim(4)], 'Color','black');

hold('on')
plot(centroidiPCA(:,1),centroidiPCA(:,2),'s','MarkerFaceColor','r','LineWidth',5)


%% Rappresentazione degli score + centroidi + (gruppi di appartenenzae etichetta provincia)
close all
idxAndRowLabels=strcat(num2str(classif),'--',etirighe);
plot(Y(:,1),Y(:,2),'o')
axislim=axis;
xlabel('Prima componente principale (Dimensione della regione)');
ylabel('Seconda componente principale (Superficie (+) e densità (-))');
text(Y(:,1),Y(:,2),idxAndRowLabels)
% Aggiungo l'asse x
line([axislim(1);axislim(2)], [0;0], 'Color','black');
% Aggiungo l'asse y
line([0;0],[axislim(3);axislim(4)], 'Color','black');

hold('on')
plot(centroidiPCA(:,1),centroidiPCA(:,2),'s','MarkerFaceColor','r','LineWidth',5)