%% ESERCIZIO I
% Caricare il dataset smoke.mat presente dentro FSDA. Costruire la tabella
% di contingenza tra le due variabili "profession" e "smoke_frequency". La
% variabile profession è qualitativa ed indica il tipo di professione. La
% variabile smoke_frequency è qualitativa e si riferisce all'attitudine
% verso il fumo. Denominare la tabella di contingenza con le prime 4
% lettere del proprio cognome (senza accenti). Effettuare un'analisi delle
% corrispondenze tra le due variabili di cui sopra e commentare il grafico.
% Aggiungere come titolo del grafico il proprio cognome ed il numero di
% matricola. Utilizzando la funzione CorAnaplot costruire il grafico di
% analisi delle corrispondenze in modo tale che il colore delle etichette
% dei punti COLONNA (attenzione punti colonna non punti riga) sia
% proporzionale alla comunalità (quota di inerzia spiegata) di ogni punto.
% Discutere la differenza in termini di inerzia spiegata tra i punti
% colonna "None" e "Medium". 

load smoke;
% Analisi delle corrispondenze partendo dalla matrice dei dati
outCA=CorAna(smoke,'datamatrix',true);
title('Riani 051485')
% Rian= tabella di contingenza tra le due variabili
Rian=outCA.Ntable;
disp('Tabella di contingenza tra le variabili "profession" e "smoke_frequency"')
disp(Rian)

% Commento al grafico di analisi delle corrispondenze
% I senior employees generalmente non fumano. I fumatori incalliti sono
% prevalentemente i Junior managers. 


% Grafico da analisi delle corrispondenze con colormap
% per i punti colonna proporzionale alla comunalità
% figure
plots=struct;
plots.ColorMapLabelCols='CntrbDim2In';
CorAnaplot(outCA,'plots',plots)

% Il punto colonna None presenta un'inerzia spiegata maggiore di Medium.
% Entrambi i punti tuttavia, hanno una comunalità vicina ad 1.
outCA.OverviewCols(["None","Medium"],["CntrbDim2In_1","CntrbDim2In_2"])

% Quota di inerzia spiegata dalle prime due dimensioni latenti di
% "None","Medium"
disp('Quota di inerzia spiegata dalle prime due dimensioni latenti dei punti colonna "None","Medium"')
disp(sum(outCA.OverviewCols{["None","Medium"],["CntrbDim2In_1","CntrbDim2In_2"]},2))


%% ESERCIZIO II
% Generare una matrice di dati di dimensione 100x4 (4 variabili X1, X2, X3
% e X4) con le seguenti caratteristiche. Correlazione tra X1 ed X2 inversa
% molto elevata (vicina a -0.9). X1 e X3 con correlazione vicino a zero.
% Correlazione diretta molto elevata tra X1 e X4. Un valore anomalo
% eccezionalmente piccolo per X4. Calcolare la matrice di correlazione per
% verificare che le 4 variabili presentino le correlazioni richieste.
% Rappresentare la matrice dei diagrammi di dispersione inserendo i boxplot
% sulla diagonale principale. Creare il grafico in coordinate parallele e
% commentarlo.
% Seed casuale per la replicabilità dei risultati
close all
rng(100)
n=100;
X1=(1:n)';
coef1=40;
% Correlazione tra X1 ed X2 inversa
% molto elevata (vicina a -0.9)
X2=2-3*X1+coef1*randn(n,1);
% correlazione tra X1 e X3 vicino a zero (X3 viene generato in maniera
% indipendente da X1)
X3=randn(n,1);
% Correlazione diretta molto elevata tra X1 e X4. Un valore anomalo
% eccezionalmente piccolo per X4
coef4=50;
X4=2+5*X1+coef4*randn(n,1);
% Aggiungo il valore anomalo
X4(1)=min(X4)-250;
Xgen=[X1 X2 X3 X4];
%Calcolare la matrice di correlazione
Rgen=corr(Xgen);
nomivars="X"+((1:4)');
Rgentable=array2table(Rgen,"RowNames",nomivars,"VariableNames",nomivars);
disp(Rgentable)

figure
spmplot(Xgen,'dispopt','box');
% Il valore anomalo univariato è evidente dal boxplot riferito a X4.

% Rappresentazione in coordinate parallele
figure
parallelplot(array2table(Xgen))
% Forte intersezione nel fascio di rette tra X1 e X2 (correlazione negativa
% molto forte)
% oppure si poteca utilizzare
%  parallelcoords(zscore(Xgen))

%% ESERCIZIO III 
% Caricare in memoria il dataset citiesItaly.mat presente in FSDA.
% Effettuare un'analisi in componenti principali utilizzando le 3 variabili
% addedval,  depos e  unemploy.
%
% addedval = Indicatore di valore aggiunto
% depos = Indicatore relativo ai depositi bancari.
% unemploy = Indicatore relativo al tasso di disoccupazione.
%
% TUTTE LE ANALISI CHE SEGUONO IN QUESTO ESERCIZIO SI RIFERISCONO AL
% DATASET CHE CONTIENE SOLO LE TRE VARIABILI  MENZIONATE SOPRA.
%
% Selezionare un numero di componenti principali che spieghino almeno il 95
% per cento della varianza totale delle tre variabili originarie.
% Interpretare le prime due componenti principali. Dall'esame del biplot,
% quali sono le province caratterizzate da elevati livelli di unemploy?
% Qual è il quadrante associato alla più elevata qualità della vita?
% Calcolare la matrice che rappresenta la miglior rappresentazione di rango
% 1 (Zhat) della matrice originaria in termini di scostamenti
% standardizzati (Z). Calcolare la somma dei quadrati delle differenze tra
% la matrice Z e la matrice Zhat. Calcolare e mostrare il contributo della
% città di Milano alla determinazione della varianza della prima componente
% latente. Calcolare e mostrare la comunalità della variabile unemploy.
% Applicare alla matrice dei dati standardizzati, il metodo di clustering
% delle k medie imponendo 2 gruppi. Utilizzare un numero di repliche pari a
% 50. Interpretare i due gruppi che sono stati ottenuti.
%
load citiesItaly.mat
Xtable=citiesItaly(:,["addedval","depos","unemploy"]);
outPCA=pcaFS(Xtable); 
disp("Tabella della quota di varianza spiegata")
disp(outPCA.explainedT)
% Per spiegare il 95% della varianza totale delle 3 variabili
% standardizzate devo prendere almeno due componenti principali

X=Xtable{:,:};
% Z= matrice degli scostamenti standardizzati
Z=zscore(X);
n=size(X,1);
% Prima comp= indicatore di povertà
% Seconda comp= valori elevati della seconda componenente principali
% indicano province con elevati depositi ma anche elevata disoccupazione.
%
% quali sono le province caratterizzate da elevati livelli di unemploy?
% Quelle che si trovano nella direzione della freccia legata ad unemploy
% (Reggio C., Enna, Catanzaro, ....)

% Calcolare la matrice che rappresenta la miglior rappresentazione di rango
% 1 (Zhat) della matrice originaria in termini di scostamenti
% standardizzati (Z). Calcolare la somma dei quadrati delle differenze tra
% la matrice Z e la matrice Zhat.

% Zhat miglior matrice di rango 1
v1=outPCA.coeff(:,1);
Zhat=Z*(v1*v1');
disp('rango Zhat')
disp(rank(Zhat))

% Somma dei quadrati delle differenze tra Z e Zhat.
sumE=sum((Z-Zhat).^2,"all");
disp('Somma dei quadrati delle differenze tra Z e Zhat')
disp(sumE)

disp('Somma dei quadrati delle differenze tra Z e Zhat')
sumEchk=sum(outPCA.explainedT{2:end,1})*(n-1);

% Verifico che sumE e sumEchk diano esattamente lo stesso risultato
assert(abs(sumE-sumEchk)<1e-12,"Attenzione la somma dei quadrati " + ...
    "delle differenze è diversa")

% Mostro il contributo assoluto della città di Milano alla determinazione della
% varianza della prima componente latente.
disp('Contributo della città di Milano alla determinazione della varianza della prima componente latente.')
disp('Score di Milano')
disp(outPCA.scoreT{"Milano",1})


%{  
    % Osservazione: dato che la varianza della prima componente latente 
    % sum(outPCA.scoreT{:,1}.^2)/(n-1)
    % è pari a lambda1 (primo autovalore)
    % il contributo relativo è uguale a 
    disp('Contributo relativo di Milano alla determinazione di lambda1')
    lambda1=outPCA.explained(1,1);
    disp(outPCA.scoreT{"Milano",1}^2/(lambda1*(n-1)))
    % Contributo dell'8 per cento circa
%}

% Comunalità di unemploy (quota di varianza della variabile unemploy
% spiegate dalle prime due componenti principali)
disp('Comunalità di Unemploy')
disp(outPCA.communalitiesT("unemploy",3))

% Cluster analysis con il metodo delle k-medie
% Si inizializza la procedura iterativa facendo
% riferimento a 50 centroidi iniziali.
[idx,Centroidi]=kmeans(Z,2,'Replicates',50);
spmplot(Xtable,idx)

% Commentare i gruppi che sono stati ottenuti. Il metodo delle k-medie
% tende a trovare gruppi di forma sferica. Un gruppo è composto da città
% caratterizzate da elevati livelli di disoccupazione, basso valore
% aggiunto e bassi depositi (città povere). L'altro gruppo è composto da
% città caratterizzate da bassi livelli  di disoccupazione, alto valore
% aggiunto e alti depositi (città ricche).


%% Domanda facoltativa
%
% Rappresentare i gruppi trovati ed i relativi centroidi nello spazio delle
% prime due componenti principali e commentare i risultati ottenuti.
% rappresentazioni dei gruppi nello spazio delle  prime due componenti
% principali
figure
Y=outPCA.score;
plot(Y(:,1),Y(:,2),'o')
% Aggiungo etichetta del nome della città e gruppo di appartenenza
text(Y(:,1),Y(:,2),strcat(Xtable.Properties.RowNames,'-',num2str(idx)))
% Aggiungo gli assi cartesiani
xline(0)
yline(0)
% Aggiungo i centroidi dei due gruppi
Cent=Centroidi*outPCA.coeff(:,1:2);
hold('on')
plot(Cent(:,1),Cent(:,2),'d','MarkerSize',30)
title('Score+indicazione dei gruppi+centroidi dei due gruppi')
xlabel('PC1')
ylabel('PC2')
% Un gruppo è caratterizzato da valori generalemnte positivi di PC1 e
% l'altro gruppo con valori prevalentemente negativi di PC1