% La zona B2:D113 del foglio dati del file Excel contiene i punteggi di tre
% test (d'intelligenza) sostenuti da 112 studenti. Importare i dati dentro
% MATLAB.
% Dopo aver standardizzato i dati:
% 1) calcolare le distanze euclidee dal centroide utilizzando le formule
% matriciali sqrt (z_i' z_i)   i=1, 2, ...., n
% 2) calcolare le distanze euclidee dal centroide utilizzando formule non
% matriciali.
% 3) Calcolare la matrice di correlazione R utilizzando le formule
% matriciali Z'Z/(n-1)
% 4) calcolare le distanze di Mahalanobis dal centroide.
% sqrt z_i R^-1 z_i , ................... i=1, 2, .... n
% Discutere se in questo esempio è preferibile lavorare con le distanze
% euclidee oppure con le distanze di Mahalanobis. Dire come sarebbero
% cambiate le distanze di Mahalanobis se avessimo operato sui dati
% originari (non standardizzati) ed avessimo utilizzato la matrice di
% covarianze anziché la matrice di correlazione. 
% (x_i -xmedio)' S^-1 (x_i-xmedio)   i=1, 2, .... n

% Ripetere  i passi 1), 2) dopo aver standardizzato i dati tramite la
% mediana ed il MAD. Sui dati standardizzati in maniera robusta calcolare
% la matrice di correlazione tramite la funzione corr. Come cambia
% la matrice di correlazione rispetto a quella calcolata al punto 3)
%
% Calcolare le
% distanze di Mahalanobis dal centroide robusto. 
% sqrt z_i R^-1 z_i ,
% ................... i=1, 2, .... n
% Costruire il boxplot delle distanze di Mahalanobis al quadrato.
% Considerare come outliers gli studenti che presentano una distanza (al
% quadrato) di Mahalanobis  superiore alla soglia  x075+3 DI dove DI è la
% differenza interquartile e x075 è il terzo quartiile delle distanze di
% Mahalanobis al quadrato.
% Rappresentare le righe della matrice dei dati tramite facce oppure
% stelle. Commentare i risultati ottenuti.
%
% Eliminare dalla matrice dei dati le righe identificate come anomale.
% Calcolare la matrice dei diagrammi di dispersione inserendo sulla
% diagonale principale i grafici ad istogrammi. Commentare i risultati
% ottenuti. 
% Sulla matrice ripulita dai valori anomali, dopo aver standardizzato i
% dati,  ridurre le dimensioni tramite la tecnica delle componenti
% principali.
% Discutere il numero appropriato di dimensioni da considerare. Calcolare
% le comunalità utilizzando solo la prima componente principale.
%
% Costuire il biplot inserendo sull'asse delle ascisse un titolo
% appropriato. Individuare lo studente più intelligente e quello meno
% intelligente nel grafico e commentare i risultati.


%% Soluzione
[X, Varnames] = xlsread('test2017B.xlsx','dati','B1:D113');

%% Standardizzo i dati
Z=zscore(X);

%% 1) calcolare le distanze euclidee dal centroide utilizzando le formule
% matriciali sqrt (z_i' z_i)   i=1, 2, ...., n
% 2) calcolare le distanze euclidee dal centroide utilizzando formule non
% matriciali.

[n,p]=size(X);
% disteuclid = vettore che contiene le distanze euclidee dal centroide di
% ogni unità
disteuclid=zeros(n,1);
for i=1:n
    disteuclid(i)=sqrt(Z(i,:)*(Z(i,:)'));
end

% disteuclidCHK = distanze euclidee dal centroide di
% ogni unità tramite formule non matriciale 
disteuclidCHK=sqrt(sum(Z.^2,2));

% Controllo che i due metodi producano gli stessi risultati
isequal(disteuclid,disteuclidCHK)

%% 3) Calcolare la matrice di correlazione R utilizzando le formule
% matriciali Z'Z/(n-1)
R=Z'*Z/(n-1);

%% 4) calcolare le distanze di Mahalanobis dal centroide.
% sqrt z_i R^-1 z_i , ................... i=1, 2, .... n
mahal=zeros(n,1);
mahalCHK=mahal;
Rinv=inv(R);
for i=1:n
mahal(i)=sqrt(Z(i,:)*Rinv*(Z(i,:)'));
% di seguito l'implementazione più efficiente suggerita da Matlab
% (FACOLTATIVO)
mahalCHK(i)=sqrt(Z(i,:)* (R\(Z(i,:)'))  );
end

% Controllo che il massimo delle differenze tra le due implementazioni sia
% un numero piccolissimo
max(abs(mahal-mahalCHK))

% distanze di Mahalanobis calcolate tramite la funzione mahalFS di FSDA 
% (FACOLTATIVO)
mahalCHK1=sqrt(mahalFS(Z,zeros(1,p),R));
max(abs(mahal-mahalCHK1))


%% Discutere se in questo esempio è preferibile lavorare con le distanze
% euclidee oppure con le distanze di Mahalanobis. Dire come sarebbero
% cambiate le distanze di Mahalanobis se avessimo operato sui dati
% originari (non standardizzati) ed avessimo utilizzato la matrice di
% covarianze anziché la matrice di correlazione. 
% (x_i -xmedio)' S^-1 (x_i-xmedio)   i=1, 2, .... n

% Risposta
% In questo esempio è preferibile lavorare con le distanze di Mahalanobis
% in quanto le variabili sono molto correlate.

% Le distanze di Mahalanobis sono invarianti per trasformazioni lineari di
% conseguenza non sarebbero cambiate se avessimo operato sui dati originari
% (non standardizzati) ed avessimo utilizzato la matrice di covarianze
% anziché la matrie di correlazione.
% Verifica empirica
% 
mu=mean(X);
S=cov(X);
mahalCHK2=sqrt(mahalFS(X,mu,S));
max(abs(mahal-mahalCHK2))


%% Ripetere  i passi 1), 2) dopo aver standardizzato i dati tramite la
% mediana ed il MAD. Sui dati standardizzati in maniera robusta calcolare
% la matrice di correlazione tramite la funzione correlazione. Come cambia
% la matrice di correlazione rispetto a quella calcolata al punto 3)
%
% Calcolare le
% distanze di Mahalanobis dal centroide robusto. 
% sqrt z_i R^-1 z_i ,
% ................... i=1, 2, .... n
% Costruire il boxplot delle distanze di Mahalanobis al quadrato.
% Considerare come outliers gli studenti che presentano una distanza (al
% quadrato) di Mahalanobis  superiore alla soglia  x075+3 DI dove DI è la
% differenza interquartile e x075 è il terzo quartiile delle distanze di
% Mahalanobis al quadrato.
% Rappresentare le righe della matrice dei dati tramite facce oppure
% stelle. Commentare i risultati ottenuti.

% Zrob = matrice che conterrà i dati standardizzati in maniera robusta
Zrob=zeros(n,p);


% Trovo i dati standardizzati tramite doppio ciclo for
% med = vettore che contiene le mediane
medians=median(X);
% mads = vettore che contiene i MAD 
mads=mad(X,1);
for i=1:n
    for j=1:p
        Zrob(i,j)=(X(i,j)-medians(j))/mads(j);
    end
end

% Metodo alternativo per trovare gli scostamenti standardizzati robusti
% tramite moltiplicazioni matriciali
Zrobchk=(X-medians)*inv(diag(mads));

% Matrice di correlazione utilizzando la funzione corr
Rchk=corr(Zrob);
% Ovviamente questa matrice è esattamente uguale alla matrice di
% correlazione R ottenuta in precedenza in quanto i coefficienti di
% correlazione lineare sono invarianti per trasformazioni lineari.
disp(max(max(abs(R-Rchk))))


%% Calcolare le
% distanze di Mahalanobis dal centroide robusto. 
% sqrt z_i R^-1 z_i ,
% ................... i=1, 2, .... n

% Il centroide robusto (ossia la mediana) dei dati standardizzati in
% maniera rousta è il vettore di 0
median(Zrob) 
% di conseguenza le distanze di Mahalanobis dal centroide robusto.
% sqrt (z_i-centrrobusto) R^-1 (z_i-centrrobusto) , diventano
% sqrt z_i R^-1 z_i ,
% ................... i=1, 2, .... n


mahalrob=sqrt(mahalFS(Zrob,zeros(1,p),R));


%% Costruire il boxplot delle distanze di Mahalanobis al quadrato.
% Considerare come outliers gli studenti che presentano una distanza (al
% quadrato) di Mahalanobis  superiore alla soglia  x075+3 DI dove DI è la
% differenza interquartile e x075 è il terzo quartiile delle distanze di
% Mahalanobis al quadrato.
% Rappresentare le righe della matrice dei dati tramite facce oppure
% stelle. Commentare i risultati ottenuti.
mahalrob2=mahalrob.^2;
boxplot(mahalrob2)

% Trovo la soglia x075+3 DI
x075=prctile(mahalrob2,75);
x025=prctile(mahalrob2,25);
soglia=x075+3* (x075-x025);
seq=1:n;
disp('Unità dichiarate outliers')
anomale=seq(mahalrob2>soglia);
disp(anomale)


%% Rappresentare le righe della matrice dei dati tramite facce oppure
% stelle. Commentare i risultati ottenuti.
glyphplot(X)
glyphplot(X,'glyph','face')
% Da questi grafici l'anomali delle unità 111 e 112 non traspare in alcun
% modo. 

%% Eliminare dalla matrice dei dati le righe identificate come anomale.
% Dato che le unità anomale si riferiscono alle ultime due righe, per
% trovare la matrice dei dati pulita è sufficiente l'istruzione che segue
Xpulita=X(1:end-2,:);

% Se le unità anomale si fossero trovate in posizione qualsiasi allora per
% trovare la matrice dei dati puliti era necessaria l'istruzione
XpulitaCHK=X;
XpulitaCHK(anomale,:)=[];
isequal(Xpulita,XpulitaCHK)



%% Calcolare la matrice dei diagrammi di dispersione inserendo sulla
% diagonale principale i grafici ad istogrammi. Commentare i risultati
% ottenuti. 
gplotmatrix(Xpulita)
% la matrice dei digrammi di dispersione mostra una distribzuione ellittica
% e una forte correlazione positiva tra le 3 variabili


%% FACOLTATIVO: scatter plot matrix utilizzando un simbolo diverso le due unità dichiarate anomale
group=ones(n,1);
group(anomale)=2;
spmplot(X,group)
% La scatter plot matrix sulle variabili orginarie rivela che le due unità
% dichiarate anomale anche se dal punto di vista univariato appartengono al
% cuore della distribuzione dal punto di vista bivariato si collocano
% sempre al margine o al di fuori della distribuzione ellittica. Queste due
% unità (studenti) sono quindi outliers multivariati!



%% Sulla matrice ripulita dai valori anomali, dopo aver standardizzato i
% dati,  ridurre le dimensioni tramite la tecnica delle componenti
% principali.
% Discutere il numero appropriato di dimensioni da considerare. Calcolare
% le comunalità utilizzando solo la prima componente principale.
%
% Costuire il biplot inserendo sull'asse delle ascisse un titolo
% appropriato. Individuare lo studente più intelligente e quello meno
% intelligente nel grafico e commentare i risultati.
Zpulita=zscore(Xpulita);
npulita=size(Zpulita,1);
% Scomposizione in valori singolari

[U,Gammastar,V]=svd(Zpulita,'econ');
% Controllo sulla svd
sqn1=sqrt(npulita-1);
Gamma=Gammastar/sqn1;
omega=1;
alpha=0;
PuntiRiga=sqn1^omega*U(:,1:2)*Gamma(1:2,1:2)^alpha;
PuntiColonna=V(:,1:2)*(Gamma(1:2,1:2)^(1-alpha))*sqn1^(1-omega);

figure
hold('on')
plot(PuntiRiga(:,1),PuntiRiga(:,2),'o')
text(PuntiRiga(:,1),PuntiRiga(:,2),num2str((1:npulita)'))

zeroes = zeros(p,1); 
quiver(zeroes,zeroes,PuntiColonna(:,1),PuntiColonna(:,2))
varlabs=Varnames;
dx=0.02;
dy=0.03;
text(PuntiColonna(:,1)+dx,PuntiColonna(:,2)+dy,varlabs,'Color','b','Interpreter','none');

axislim=axis;
% 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');

% Commenti:
% La prima componente è determinata da tutti e tre i test con pesi
% pressoché uguali.
disp(V(:,1))
% Dato che i segni sono tutti negativi, la prima componente si interpreta
% come un indicatore latente di mancanza di intelligenza
xlabel('Mancanza di intelligenza')

% Quota di varianza spiegata dalla prima CP
disp(Gamma(1,1)^2/p)

% Quota di varianza spiegata dalle prime due CP
disp((Gamma(1,1)^2+Gamma(2,2)^2)/p)

% criterio soglia =.95^p=
disp(['Criterio soglia=' num2str(0.95^p)])

% Utilizzando il criterio soglia =.95^p ,dato che la prima componente
% principale spiega una quota molto superiore del valore soglia ha senso
% considerare solo la prima componente principale.

% Il biplot viene costruito solo a fini didattici in quanto per
% sintetizzare le informazioni è sufficiente la prima componente principale

% L'angolo molto vicino a 90 gradi tra la seconda componente e le tre
% frecce ribadisce che la seconda componente è poco correlata con le tre
% variabili in esame

% Calcolare
% le comunalità utilizzando solo la prima componente principale.
% Analisi delle comunalità 
disp((V(:,1)*Gamma(1,1)).^2)
% La quota di varianza di ciascuno dei tre test spiegata dalla prima
% componente è al di sopra del 93.5%

% Lo studente più intelligente è quello che ha lo score più basso della
% prima componente principale. Studente numero  9
% Lo studente meno intelligente è quello che ha lo score più alto della
% prima componente principale. Studente numero 35