%% Esercizio
% Reinizializzare a default il generatore dei numeri casuali
% Generare una matrice di dimension n1-by-p (n1=10 e p=2) di numeri casuali
% dalla distribuzione N(0,1). Generare una matrice di dimension n2-times-p
% (n2=20 e p=2) di numeri casuali dalla distribuzione N(0+dx,1) con dx=5.
% Mostrare i due gruppi di unità (utilizzando simboli diversi per i due
% gruppi) in un diagramma a dispersione.
% Combinare i due gruppi in una matrice X di dimensione (n1+n2)-by-p;
% Costruire un grafico a dispersione delle due colonne della matrice X con
% l'aggiunta del numero di ogni unità
%
% CLUSTERING GERARCHICO
%
% Applicare la metodologia del clustering gerarchico (utilizzando il
% metodo del legame medio e la metrica della distanza euclidea). Costruire
% il relativo dendrogramma, mostrando tutti i passi della procedura di
% agglomerazione.
% Tagliare il dendrogramma alla soglia di distanza 3.5.
% Personalizzare l'orientamento del dendrogramma  e/o la 'ColorThreshold'.
% Opzioni 'Orientation', e 'ColorThreshold',
% 
% CLUSTERING NON GERARCHICO
%
% Applicare la metodologia del clustering non gerarchico utilizzando la
% metrica della distanza della città a blocchi e il metodo di clustering
% delle kmedie utilizzando 3 gruppi. Rappresentare la classficazione
% tramite la funzione spmplot
%
% Esplorare come cambia la classificazione utilizzando un numero di
% gruppi pari a k (con k=2, 3, 4 e 5). Visualizzare il risultato della
% classificazione tramite la funzione gscatter. Utilizzare un grafico a 4
% pannelli. Il pannello in alto a sinistra contiene la classificazione con
% 2 gruppi, il pannello in alto a destra contiene la classificazione con 3
% gruppi ....
% Ripetere la procedura aggiungendo ai 4 pannelli i centroidi dei gruppi.
%
%
% Tramite la funzione di MATLAB evalclusters scegliere il numero ottimo di gruppi.
% Come numero di gruppi da esplorare scegliere k=1:6
% Come criterio da utilizzare scegliere CalinskiHarabasz
% Osservazione: criterio di CalinskiHarabasz= devianza tra i gruppi / devianza nei gruppi
% Rappresentare graficamente l'indice di CalinskiHarabasz in funzione di k
% Commentare il grafico.

% ESEMPI PRATICI UTILIZZANDO DATASETS NOTI NELLA LETTERATURA DELLA CLUSTER
% ANALYSIS

% Fare il clear della memoria
% Caricare in memoria il dataset contenuto nel file geyser.txt dell'FSDA
% toolbox.
% Per ulteriori informazioni su questo dataset di veda la sezione datasets
% dell'FSDA toolbox.
% Costruire il grafico a dispersione delle due
% variabili nel dataset. Dall'analisi grafica quanti gruppi sembrano essere
% presenti? Ci sono outliers?
% Esplorare i risultati che si ottengono utilizzando la metodologia delle
% k-medie e la metrica della città a blocchi k=3
% Commentare i risultati ottenuti e l'effetto degli outliers sulla
% classificazione. Visualizzare i risultati della classificazione
% utilizzando la funzione spmplot.
%
% Esplorare come varia la classificazione utilizzando diverse metriche.
% Fissare k =3 ed esplorare (in un grafico a 4 pannelli) la classificazioen
% che si ottiene con le seguenti 4 metriche
% 'sqeuclidean', 'cityblock', 'cosine', 'correlation'
%

%
% Calcolare il numero ottimo dei gruppi che viene suggerito
% dall'applicazione della funzione evalclusters (utilizzare la distanza Euclidea).
%

% Utilizzando la funzione tkmeans (trimmed k-means) dell'FSDA toolbox,
% analizzare la classificazione che si ottiene quando si impone un numero
% di gruppi k=3 e una frazione di trimming alpha pari a 0.05.
% Analizzare la stabilità dei risultati al variare della frazione di
% trimming alpha.

% Osservazione: finora abbiamo lavorato con gruppi di forma sferica.
% L'obiettivo dell'esempio che segue è quello di mostrare la necessità di
% lavorare con le distanze di Mahalanobis al posto delle distanze Euclidee.
%
%
% Resettare il generatore di numeri casuali.
% Generare un numero di unità n1 (n1=100) da una distribuzione normale
% bivariata con vettore delle medie mu=(2 3)' e matrice di covarianze
% 1    1.5
% 1.5   3
% Generare un numero di unità n2 (n2=100) da una distribuzione normale
% bivariata con vettore delle medie mu=(5 3)' e matrice di covarianze
% 1    1.4
% 1.4   3
% Combinare i due gruppi in una matrice di dimensione n1+n2-by-2.
% Costruire il diagramma di dispersione dei due gruppi che sono stati generati

% Applicare la metodologia delle k-medie utilizzando la metrica della
% distanza Euclidea e un numero di gruppi pari a 2. Commentare i risultati
% della classificazione che si ottiene
% Calcolare il numero ottimo di gruppi che risulta dall'applicazione della funzione
% evalclusters (come numero minimo e massimo di k scegliere 1 e 10).
% Analizzare la classificazione che si ottiene (tramite il metodo kmeans e
% la metrica della distanza Euclidea) inserendo il numero ottimo di gruppi
% trovato in automatico dalla funzione evalclusters.
% Allo scopo di ottenere risultati stabili, utilizzare un numero di
% repliche pari a 100 quando si chiama la funzione kmeans.
% Applicare al dataset la funzione tclust utilizzando un numero di gruppi
% pari a 2, una percentuale di trimming pari a zero e un restriction factor uguale a 100.
% Visualizzare i risultati della classificazione e gli ellissi di
% confidenza associati a ciascun gruppo.
% Commentare i risultati ottenuti.

% Fare il clear della memoria e chiudere tutti i grafici
% Caricare in memoria i dati della qualità della vita
% Applicare il metodo di raggruppamento non gerarchico tclust utilizzando
% un numero di gruppi pari a 3, una frazione di trimming pari al 3 per
% cento ed un valore del fattore di restrizione pari a 10.
% Per avere la replicabilità dei risultati impostare l'opzione nsamp a 1000
% Visualizzare e commentare i risultati della classificazione
% Mostrare le unità (province che non sono state classificate in quanto
% anomale)


%% Generazione di due gruppi
% Reinizializzare a default il generatore dei numeri casuali
% Generare una matrice di dimension n1-by-p (n1=10 e p=2) di numeri casuali
% dalla distribuzione N(0,1). Generare una matrice di dimension n2-times-p
% (n2=20 e p=2) di numeri casuali dalla distribuzione N(0+dx,1) con dx=5.
% Mostrare i due gruppi di unità (utilizzando simboli diversi per i due
% gruppi) in un diagramma a dispersione.
% Combinare i due gruppi in una matrice di dimensione (n1+n2)-by-p;

rng default
n1=10;
n2=20;
p=2;
% dx = amount of displacement
dx=5;
X1=randn(n1,p);
X2=dx+randn(n2,p);
% Mostrare i due gruppi di unità (utilizzando simboli diversi per i due
% gruppi) in un diagramma a dispersione.
hold('on')
plot(X1(:,1),X1(:,2),'bo')
plot(X2(:,1),X2(:,2),'rs')

% Combinare i due gruppi in una matrice di dimensione (n1+n2)-by-p;
X=[X1; X2];

% Costruire un grafico a dispersione delle due colonne della matrice X con
% l'aggiunta del numero di ogni unità
figure
lab=num2str((1:n1+n2)');
plot(X(:,1),X(:,2),'o')
text(X(:,1),X(:,2),lab)

%% Cluster analysis gerarchica
% Applicare la metodologia del clustering gerarchico (utilizzando il
% metodo del legame medio e la metrica della distanza euclidea). Costruire
% il relativo dendrogramma, mostrando tutti i passi della procedura di
% agglomerazione.
% Personalizzare l'orientamento del dendrogramma  e/o la 'ColorThres
% Metodo gerarchico del legamo medio.
close all

% Es. metodo del legame medio ('average')
% Metrica utilizzata: distanza euclidea
tree = linkage(X,'average','euclidean');

% Creazione del dendrogramma
% Il secondo argomento pari a 0
% specifica di mostrare tutti i nodi
dendrogram(tree,0)
% dendrogram(tree,0,'ColorThreshold',3.5)

% Taglio del dendrogramma  alla soglia di distanza 3.5
group = cluster(tree,'cutoff',3.5,'Criterion','distance');

% spmplot di FSDA per visualizzare i gruppi ottenuti
spmplot(X,group);

%% Dendrogrammi personalizzati
% Personalizzare l'orientamento del dendrogramma  e/o la 'ColorThreshold'.
tree = linkage(X,'average');
H = dendrogram(tree,'Orientation','left','ColorThreshold','default');
% Impostazione della linewidth del dendrogramma
set(H,'LineWidth',2)


%% Metodo non gerarchici
% CLUSTERING NON GERARCHICO
%
% Applicare la metodologia del clustering non gerarchico utilizzando la
% metrica della distanza della città a blocchi e il metodo di clustering
% delle kmedie utilizzando 3 gruppi. Rappresentare la classficazione
% tramite la funzione spmplot

% Ad esempio metodo delle k medie
k=3;
idx = kmeans(X,k,'Distance','cityblock');
spmplot(X,idx)


%% Analisi di come varia la classificazione al variare di k
% Esplorare come cambia la classificazione utilizzando un numero di
% gruppi pari a k (con k=2, 3, 4 e 5). Visualizzare il risultato della
% classificazione tramite la funzione gscatter. Utilizzare un grafico a 4
% pannelli. Il pannello in alto a sinistra contiene la classificazione con
% 2 gruppi, il pannello in alto a destra contiene la classificazione con 3
% gruppi ....
for k=2:5
    subplot(2,2,k-1)
    idx = kmeans(X,k,'Distance','cityblock');
    gscatter(X(:,1),X(:,2),idx,[],[],25)
    title(['k=' num2str(k)])
end

%% Analisi di come varia la classificazione al variare di k (con aggiunta dei centroidi)
%
% Ripetere la procedura aggiungendo ai 4 pannelli i centroidi dei gruppi.
%
for k=2:5
    subplot(2,2,k-1)
    idx = kmeans(X,k,'Distance','cityblock');
    gscatter(X(:,1),X(:,2),idx,[],[],25)
    Centroidi=grpstats(X,idx);
    hold('on')
    gscatter(Centroidi(:,1),Centroidi(:,2),1:k,[],'p',35)
    title(['k=' num2str(k) ' con i centroidi dei gruppi'])
end


%% Scelta del numero ottimo di gruppi
% Tramite la funzione di MATLAB evalclusters scegliere il numero ottimo di gruppi.
% Come numero di gruppi da esplorare scegliere k=1:6
% Come criterio da utilizzare scegliere CalinskiHarabasz
% Osservazione: criterio di CalinskiHarabasz= devianza tra i gruppi / devianza nei gruppi
% Rappresentare graficamente l'indice di CalinskiHarabasz in funzione di k
% Commentare il grafico.
figure
eva = evalclusters(X,'kmeans','CalinskiHarabasz','KList',1:6);
plot(eva);
% Questo grafico mostra chiaramente che il numero ottimo dei gruppi è pari
% a 2 

%% Geyser data
close all
Y=load('geyser2.txt');
plot(Y(:,1),Y(:,2),'o')

% Commento: il grafico rivela la presenza di 3 gruppi (sferici) e di alcune unità
% anomale concentrate nella zona in basso a sinistra del diagramma di
% dispersione.

%% Metodo delle k medie con diverso numero di gruppi
% Esplorare la classificazione che si ottiene utilizzando rispettivamente
% k=2, k=3 e k=4
k=3;
% idx = kmeans(Y,k,'Distance','sqeuclidean');
idx = kmeans(Y,k,'Distance','cityblock');
spmplot(Y,idx)

% Commento: la presenza di outliers ha un effetto enorme sulla
% classificazione.

%% Analisi della stabilità dei risultati al variare della metrica
tipodistanza={'sqeuclidean', 'cityblock', 'cosine', 'correlation'};
k=3; % k viene sempre fissato a 3
close all
for j=1:4
    subplot(2,2,j)
    idx = kmeans(Y,k,'Distance',tipodistanza{j});
    gscatter(Y(:,1),Y(:,2),idx,[],[],25)
    Centroidi=grpstats(Y,idx);
    hold('on')
    h=gscatter(Centroidi(:,1),Centroidi(:,2),1:k,[],'*',35);
    
    title(['k=' num2str(k) ' distanza' tipodistanza{j}])
end


%% Choose the best number of clusters (in presenza di un dataset contenente outliers)
close all
eva = evalclusters(Y,'kmeans','CalinskiHarabasz','KList',1:6);
plot(eva)
% Il numero ottimo di gruppi che viene suggerito è 4.
% La classificazione con 4 gruppi risulta inaccettabile in quanto
% fortemente influenzata dalla presenza di valori anomali

%% Trimmed k-means
k=3;
% Nella metodologia del tkmeans si deve specificare la frazione alpha di
% osservazioni che non verranno classificate
% (ossia che verranno allocate al gruppo 0).
alpha=0.10;
out=tkmeans(Y,k,alpha,'plots',1);



%% Generazione di due gruppi ellittici
% Resettare il generatore di numeri casuali.
% Generare un numero di unità n1 (n1=100) da una distribuzione normale
% bivariata con vettore delle medie mu=(2 3)' e matrice di covarianze
% 1    1.5
% 1.5   3
% Generare un numero di unità n2 (n2=100) da una distribuzione normale
% bivariata con vettore delle medie mu=(5 3)' e matrice di covarianze
% 1    1.4
% 1.4   3
% Combinare i due gruppi in una matrice di dimensione n1+n2-by-2.
rng default  
n1=100;
mu = [2,3];
cov1=1.5;
sigma = [1,cov1;cov1,3];
X1 = mvnrnd(mu,sigma,n1);

mu = [5,3];
cov2=1.4;
n2=100;
sigma = [1,cov1;cov1,3];
X2 = mvnrnd(mu,sigma,n2);
X=[X1; X2];

%% Diagramma di dispersione dei due gruppi che sono stati generati
close all
hold('on')
plot(X1(:,1),X1(:,2),'bo')
plot(X2(:,1),X2(:,2),'rs')

%% Metodologia delle k medie
% Applicare la metodologia delle k-medie utilizzando la metrica della
% distanza Euclidea e un numero di gruppi pari a 2. Commentare i risultati
% della classificazione che si ottiene
% Calcolare il numero ottimo di gruppi che risulta dall'applicazione della funzione
% evalclusters.

idx = kmeans(X,2,'Distance','sqeuclidean');
spmplot(X,idx)

% Commento: la metodologia delle k medie permette di trovare solo gruppi di
% forma sferica ma non gruppi di forma ellissoidale.

%% Funzione evalcluster in presenza di gruppi ellittici
close all
% Calcolare il numero ottimo di gruppi che risulta dall'applicazione della funzione
% evalclusters (come numero minimo e massimo di k scegliere 1 e 10).

eva = evalclusters(X,'kmeans','CalinskiHarabasz','KList',1:10);
plot(eva)

% Il numero ottimo di gruppi suggerito è pari a 8.
% Il numero vero di gruppi era pari a 2!!!!

%% Analisi classificazione con k ottimo secondo evalclusters
% Analizzare la classificazione che si ottiene (tramite il metodo kmeans e
% la metrica della distanza Euclidea) inserendo il numero ottimo di gruppi
% trovato in automatico dalla funzione evalclusters.
% Allo scopo di ottenere risultati stabili, utilizzare un numero di
% repliche pari a 100 quando si chiama la funzione kmeans.

idx = kmeans(X,8,'Distance','sqeuclidean','Replicates', 100);
spmplot(X,idx)
% Commento: i gruppi che vengono trovati hanno una forma sferica.
% Il primo gruppo viene suddiviso in 4 sottogruppi. La stessa cosa accade
% per il secondo.


%% tclust
% Applicare al dataset la funzione tclust utilizzando un numero di gruppi
% pari a 2, una percentuale di trimming pari a zero e un restriction factor uguale a 100.

out = tclust(X,2,0,10,'plots',1);
% Il quarto input della funzione tclust controlla il massimo rapporto
% ammesso tra la lunghezza del più grande asse maggiore di qualsiasi ellisse associato a
% ciascun gruppo e il più piccolo asse minore di qualsiasi ellisse
% associato a ciascun gruppo. restrfactor=1 significa che vogliano trovare
% solo gruppi di forma sferica. Osservazione: nella sfera il rapporto tra
% la lunghezza degli assi è pari a 1.
% Un altro argomento opzionale è equalweights: nella funzione obiettivo
% vogliamo o meno pesare le numerosità dei gruppi (ossia dare un peso
% maggiore ai gruppi più numerosi). Nella metodolgia delle k medie questo
% peso non si assegna (cioè equalweights=1). Dentro tclust il default di
% equalweights è 0 (ossia i gruppi più numerosi pesano di più).

% Ad esempio la chiamata con
% 1) livello di trimming pari a zero (terzo argomento di input)
% 2) restriction factor uguale a 1 (quarto argomento di input)
% 3) 'name',value 'equalweights',1
% equivale ad utilizzare il metodo delle
% k-medie e si trovano solo gruppi di forma sferica.
% out = tclust(X,2,0,1,'plots',1,'equalweights',1);
% spmplot(X,out.idx)

%% Visualizzazione della classificazione aggiungendo ellissi di confidenza 
% Visualizzare i risultati della classificazione e gli ellissi di
% confidenza associati a ciascun gruppo.
% Commentare i risultati ottenuti.
overl=struct;
overl.type='ellipse';
spmplot(X,'group',out.idx,'overlay',overl)

% Commento: la funzione tclust è basata sulle distanze di Mahalanobis e
% consente di ottenere anche gruppi di forma ellittica.


%%  Clustering robusto dei dati della qualità della vita

% Fare il clear della memoria e chiudere tutti i grafici
clear 
close all


% Caricare in memoria i dati della qualità della vita
Xtable=readtable('benessere.xlsx','Sheet','X (database originale)','ReadRowNames',true,'Range','A3:H106');
% X = matrice di double senza nomi delle righe e nomi delle colonne
X=table2array(Xtable);

% varalbs = cell che contiene i nomi delle varaibili
varlabs=Xtable.Properties.VariableNames;
nomiprov=Xtable.Properties.RowNames;
[n,p]=size(X);
Z=zscore(X);


% Applicare il metodo di raggruppamento non gerarchico tclust utilizzando
% un numero di gruppi pari a 3, una frazione di trimming pari al 3 per
% cento ed un valore del fattore di restrizione pari a 10.
% Per avere la replicabilità dei risultati impostare l'opzione nsamp a 1000

out=tclust(Z,3,0.03,10,'plots',0,'nsamp',1000);

%% Visualizzare e commentare i risultati della classificazione
close all
spmplot(Xtable,out.idx)

% Mostrare le unità (province che non sono state classificate in quanto
% anomale)
% Unità non classificate in quanto valori anomali (trimmate)
disp('Province non classificate in quanto anomale')
disp(nomiprov(out.idx==0))

% gscatter(X(:,1),X(:,2),out.idx)