%% TESTO
% 1) Importare in formato table (utilizzare come nomi delle righe le
% province) i dati presenti nel file benessere.xlsx
% 2) Calcolare la matrice di correlazione, la matrice che contiene le
% componenti principali (Y=Z V). Qual è la matrice di covarianze di Y? 
% 3) Calcolare la quota di varianza (relativa e cumulata) spiegata da ogni
% componente principale. Visualizzare la varianza spiegata dalle diverse
% componenti tramite diagramma di Pareto
% 4) Calcolare le correlazioni tra le variabili originarie e le
% componenti principali (matrici di componenti).  Calcolare le correlazioni
% direttamente applicando la formula $V \sqrt \Lambda$ oppure facendo
% riferimento alla funzione corr. 
% Verificare che la somma dei quadrati delle colonne della matrice di
% componenti sono uguali ai rispettivi autovalori. 
% Rappresetare tramite grafici a barre la correlazione tra le variabili
% originarie e le due componenti principali. Inserire i due grafici in due
% subplots. Aggiungere nel grafico a barre sopra le barre l'indicazione del
% valore numerico.
% 5) Calcolare le componenti principali (non standardizzate) Y=Z*V;
% Verificare tramite la funzione corr che le correlazioni tra la j-esima
% componente principale e le p variabili originarie sono uguali alla
% colonna della j-esima della matrice di componenti. Verificare che la
% matrice di covarianze delle p componenti principali è diagonale e
% contiene sulla diagonale principale gli autovalori di R. 
% 6) Partendo dalla matrice di componenti, calcolare le  comunalità
% (quote di varianza di ogni variabile spiegate dalla prime due componenti
% principali)
% Verificare che le comunalità (indici R2) si possono anche ottenere
% tramite una regressione multipla in cui la variabile in oggetto
% (variabile dipendente) viene regredita sulle prime due componenti
% principali (variabili esplicative). Verificare questo risultato con
% riferimento alla seconda variabile.
% 7) Rappresentare nel piano cartesiano tramite frecce delle prime due colonne
% della matrice di componenti, prima utilizzando la funzione biplot, poi
% tramite l'implementazione manuale e successivamente tramite la funzione
% quiver
% 8) Rappresentare nel piano cartesiano gli score non standardizzati (prime
% due colonne della matrice Z*V inserendo le etichette di riga) e poi gli score
% standardizzati (componenti principali standardizzate = componenti
% principali con matrice di covarianze uguale all'identità), inserendo le
% etichette appropriate per gli assi cartesiani.
% Verificare che la matrice degli score standardizzati si può ottenere tramite
% la scomposizione in valori singolari della matrice Z passando per la matrice U.


%% Soluzione

%%  CARICAMENTO DATI
% 1) Importare in formato table (utilizzare come nomi delle righe le
% province) i dati presenti nel file benessere.xlsx
Xtable=readtable('benessere.xlsx','ReadRowNames',true);
% X = matrice di double senza nomi delle righe e nomi delle colonne
X=table2array(Xtable);

% nameXvars = cell che contiene i nomi delle variabili
nameXvars=Xtable.Properties.VariableNames;

[n,p]=size(X);

%% 2) Calcolare la matrice di correlazione, 
% Standardizzo i dati
Z=zscore(X);
% calcolo matrice di correlazione
R=cov(Z);
% mostro la matrice di correlazione in formato table
Rtable=array2table(R,'RowNames',nameXvars,'VariableNames',nameXvars);
format bank
disp(Rtable)
format short
% R poteva essere ottenuta direttamente da X utilizzando la funzione corr

%% Autovettori e autovalori della matrice di correlazione
[V,La]=eig(R);
la=diag(La);
[aa,indsor]=sort(la,'descend');
% Riordino le colonne della matrice V e La in modo tale che
% La(1,1) sia il grande autovalore e V(:,1) sia l'autovettore associato
% La(2,1) sia il secondo più grande autovalore e V(:,2) sia l'autovettore associato
V=V(:,indsor);
lasor=la(indsor);
La=diag(lasor);

% Controllo la traccia di La sia uguale alla traccia di R
disp('Traccia della matrice diagonale degli autovalori')
disp(trace(La))
la=diag(La);

%% CALCOLO COMPONENTI PRINCIPALI
% Calcolare la matrice che contiene le
% componenti principali (Y=Z V) 
Y=Z*V;

% La matrice di varianze e covarianze di Y è la matrice \Lambda
% \Lambda = \Gamma * \Gamma
% matrice degli autovalori della matrice di correlazione
% Le componenti principali sono tra loro incorrelate
% In altri termini la corr(Y) mi produce la matrice identità
disp('Matrice di correlazione calcolata sulle CP')
corr(Y)

% Calcolo della matrice di covarianze di Y
% covY contiene sulla diagonale principale gli autovalori
covY=cov(Y);

% Domanda: qual è la matrice di covarianze di Y? 
% Risposta:le componenti principali sono tra loro incorrelate di
% conseguenza la matrice di covarianze di Y è diagonale.
% Sulla diagonale principale ci sono gli autovalori.

namePCs=cellstr([repmat('PC',p,1) num2str((1:p)')]);
covYtable=array2table(covY,'RowNames',namePCs,'VariableNames',namePCs);
disp('Matrice di varianze e covarianze delle CP')
disp(covYtable)

% Controllo che la somma di diag(covY) = p
disp('Somma degli autoalori=p=traccia(R)')
disp(sum(diag(covY)))

%% VARIANZA SPIEGATA
% 3) Calcolare la quota di varianza (relativa e cumulata) spiegata da ogni
% componente principale
% Calcolo la quota di varianza spiegata da ciascuna CP
autoval=[lasor 100*(lasor)/p 100*cumsum(lasor)/p];
namecols={'Autovalori' 'Var_spiegata' 'Var_cum_spiegata'};
autovaltable=array2table(autoval,'RowNames',namePCs,'VariableNames',namecols);


%% VISUALIZZAZIONE VARIANZA SPIEGATA TRAMITE DIAGRAMMA DI PARETO
figure()
pareto(autoval(:,1),namePCs)
xlabel('Componenti principali')
ylabel('Varianza spiegata (%)')


%% INTERPRETAZIONE DELLE COMPONENTI
% 4) Calcolare le correlazioni tra le variabili originarie e le prime due
% componenti principali (matrici di componenti). 
MatrComp=V*sqrt(La);



%% Ottengo direttamente MatrComp tramite la funzione corr 
% e verifico che il risultato sia lo stesso
MatrCompchk=corr(Z,Y);
assert(max(abs(MatrComp-MatrCompchk),[],'all')<1e-10,"Errore di programmazione" + ...
    "Le correlazioni tra le variabili originarie e le PC non sono corrette")


% MatrComp ha dimensione $p \times p$. Nel nostro esempio 7*7;
% Ora la ridefinisco prendendo solo le correlazioni tra le variabili
% originarie e le prime due componenti principali
MatrComp=MatrComp(:,1:2);

%% Matrice di componenti in formato table

MatrCompt=array2table(MatrComp,"RowNames",nameXvars,"VariableNames",namePCs(1:2));
disp('Correlazioni tra le variabili originarie e le CP')
disp(MatrCompt)

%% SOMMA QUADRATI COLONNE DELLA MATRICI DI COMPONENTI
% Verificare che la somma
% dei quadrati delle colonne della matrice di componenti sono uguali ai
% rispettivi autovalori.

% La somma dei quadrati della colonna j della matrice di componenti è
% pari al j-esimo autovalore
j=1;
disp(['La somma dei quadrati della colonna ' num2str(j) ' della matrice di componenti'])
sum(MatrComp(:,j).^2)
disp(['è uguale all'' autovalore ' num2str(j) '=' num2str(lasor(j))])


%% Visualizzazione grafica matrici di componenti
% Rappresetare tramite grafici a barra la correlazione tra le variabili
% originarie e le due componenti principali. Inserire i due grafici in due
% subplots. Aggiungere nel grafico a barra sopra le barre l'indicazione del
% valore numerico.
%
% i= componente da visualizzare nel pannello in alto
i=1; 
% j= componente da visualizzare nel pannello in basso
j=2; 
% Le label sull'asse x del grfico a barre sono i nomi delle variabili, il
% secondo argomento di catergorical specifica che vogliano usare l'ordine
% delle labels speecificato dentro varlabs e non l'ordine alfabetico
xlabels=categorical(nameXvars,nameXvars);
figure
subplot(2,1,1)
b=bar(xlabels, MatrComp(:,i),'g');
title(['Correlazioni con PC' num2str(i)])

% Aggiungere nel grafico a barra sopra le barre l'indicazione del
% valore numerico.
% xtips e ytips sono le coordinate numeriche dove inserire le etichette del
% valore numerico
xtips = b.XEndPoints;
ytips = b.YEndPoints;
barlabels = string(round(MatrComp(:,1),2));
text(xtips,ytips,barlabels,'HorizontalAlignment','center',...
    'VerticalAlignment','bottom')

subplot(2,1,2)
b=bar(xlabels, MatrComp(:,j),'g');
title(['Correlazioni con PC' num2str(i)])

% Aggiungere nel grafico a barra sopra le barre l'indicazione del
% valore numerico.
xtips = b.XEndPoints;
ytips = b.YEndPoints;
barlabels = string(round(MatrComp(:,j),2));
text(xtips,ytips,barlabels,'HorizontalAlignment','center',...
    'VerticalAlignment','bottom')
title(['Correlazioni con PC' num2str(j)])


%% 5) Verificare tramite la funzione corr le correlazioni tra la prima
% componente princiale e le p variabili originarie
% Y = matrice che contiene le componenti principali
% Prima colonna di Y = prima componente principale
% Seconda colonna di Y = seconda componente principale
% ....

j=1;
MatrcompCHK=corr(Y(:,j),Z)';
disp(['Correlazioni tra la CP' num2str(j) 'e le variabili originarie'])
disp('Differenza tra le due implementazioni')
disp(max(abs(MatrComp(:,j)-MatrcompCHK)))


%% COMUNALITA'
% 6) Partendo dalla matrice di componenti, calcolare le  comunalità
% (quote di varianza di ogni variabile spiegate dalla prime due componenti
% principali)
disp('Comunalità: quote di varianza di ogni variabile spiegate dalle prime due CP')
Comu=sum(MatrComp.^2,2);

Comutable=array2table(Comu,'RowNames',nameXvars,...
    'VariableNames',{'Communalities'});
disp(Comutable)


%% VERIFICA COMUNALITA' TRAMITE REGRESSIONE LINEARE
% Verificare che le comunalità (indici R2) si possono anche ottenere
% tramite una regressione multipla in cui la variabile in oggetto
% (variabile dipendente) viene regredita sulle prime due componenti
% principali (variabili esplicative). Verificare questo risultato con
% riferimento alla seconda variabile ("depositi").
out=fitlm(Y(:,1:2),X(:,2));
disp('R2 calcolato tramite la regressione lineare')
disp(out.Rsquared.Ordinary)
disp('Comunalità della seconda variabile')
disp(Comutable(2,1));


%% 7) Rappresentare nel piano cartesiano tramite frecce le prime due colonne
% della matrice di componenti utilizzando la funzione biplot e poi tramite
% l'implementazione manuale

%% Rappresentazione grafica delle correlazioni tramite la funzione biplot
biplot(MatrComp(:,1:2),'varlabels',nameXvars)
%   Si noti che
%   BIPLOT imposes a sign convention, forcing the element with largest
%   magnitude in each column of COEFS to be positive.  This flips some of the
%   vectors in COEFS to the opposite direction, but often makes the plot
%   easier to read.  Interpretation of the plot is unaffected, because
%   changing the sign of a coefficient vector does not change its meaning.

%% Rappresentazione grafica correlazioni tramite implementazione manuale 
close all
zeroes = zeros(p,1); 
arx = [zeroes MatrComp(:,1) ]';
ary = [zeroes MatrComp(:,2)]';
% Aggiungi le linee che si dipartono dall'origine
line(arx,ary, 'Color','b', 'LineStyle','-',  'Marker','none')
% Aggiuni il simbolo di marcatore (il puntino) alla fine di ogni linea 
line(arx(2,:),ary(2,:), 'Color','k', 'Marker','.', 'LineStyle','none')

dx=0.02;
dy=0.03;
text(MatrComp(:,1)+dx,MatrComp(:,2)+dy,nameXvars);
% Aggiunge la griglia
grid on;
% Plot axes and label the figure.
axlimHi = 1.1*max(abs(MatrComp(:)));
axlimLo = -axlimHi;
% Istruzione per aggiungere gli assi con una riga sola
line([axlimLo axlimHi NaN 0 0],[0 0 NaN axlimLo axlimHi], 'Color','black');
xlabel('Prima componente principale (Indice di povertà)');
ylabel('Seconda componente principale (Indice di malessere delle aziende)');
% axis tight= Fit the axes box tightly around the data by setting the
% axis limits equal to the range of the data.
axis tight
% Se si vuole la stessa unità di misura per i due assi
% axis equal

%% Rappresentazione delle correlazioni utilizzando la funzione quiver
% In alternativa le frecce si possono tracciare utilizzando la funzione
% quiver
quiver(arx(1,:)',ary(1,:)',arx(2,:)',ary(2,:)')
hold('on')
% viene aggiunto il testo
text(MatrComp(:,1)+dx,MatrComp(:,2)+dy,nameXvars);

% Rimuovere il commento che segue se si desidera far passare gli assi per
% l'origine
% ax = gca;
% ax.XAxisLocation = 'origin';
% ax.YAxisLocation = 'origin';

% Per aggiungere la linea con coordinate x -1 1
% e coordinate y 0 0 (asse x) l'istruzione è 
% line([-1;1],[0;0])

% Istruzione per aggiungere gli assi con una riga sola
line([axlimLo axlimHi NaN 0 0],[0 0 NaN axlimLo axlimHi], 'Color','black');
xlabel('Prima componente principale (Indice di povertà)');
ylabel('Seconda componente principale (Indice di malessere delle aziende)');
% axis tight= Fit the axes box tightly around the data by setting the
% axis limits equal to the range of the data.
axis tight
% Se si vuole la stessa unità di misura per i due assi
% axis equal


%% Punto 8) 
% Rappresentare nel piano cartesiano gli score non standardizzati (prime
% due colonne della matrice inserendo le etichette di riga) e poi gli score
% standardizzati (componenti principali standardizzate = componenti
% principali con matrice di covarianze uguale all'identità), inserendo le
% etichette appropriate per gli assi cartesiani.
% Verificare che la matrice degli score standardizzati si può ottenere tramite
% la scomposizione in valori singolari della matrice Z passando per la matrice U.

%% Rappresentazione nel piano cartesiano degli scores (non standardizzati)
% chiudo preliminarmente tutti i grafici finora calcolati
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),Xtable.Properties.RowNames)
% 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');

%% Calcolo delle componenti principali standardizzate
% Yst= componenti principali standardizzate cov(Yst) = matrice identità
Yst=zscore(Y);

% In maniera matrice Yst può essere calcolata come
% Latex notation %Y \times \Gamma^{-1}%
Ystchk=Y*sqrt(inv(La));

% Controllo che gli score standardizzati si possono ottenere direttamente
% tramite la svd di Z (a meno eventualmente del segno)
% Verificare la matrice degli score standardizzati si può ottenere tramite
% la scomposizione in valori singolari della matrice Z.
[U,Gamma,Vchk]=svd(Z,'econ');
YstCHK=sqrt(n-1)*U;
Gamma=Gamma/sqrt(n-1);

% \Gamma*\Gamma = matrice degli autovalori della matrice di correlazione
Lachk=Gamma*Gamma;


% La matrice Vchk coincide con la matrice V trovata in precedenza
% attraverso la funzione eig. Si notiche che i segni di qualche colonna
% della matrice Vchk potrebbero essere scambiati rispetto a V
% In altri termini: può capitare che  V(:,j)=-Vchk(:,j)
%(v. ad esempio colonne 4:7)


% Differenza tra i due metodi di ottenere gli score standardizzati
% Prendo per entrambe le matrici la funzione abs poiché i segni potrebbero
% essere scambiati 
disp(max(max(abs(abs(Yst)-abs(YstCHK)))))

%% Rappresentazione nel piano cartesiano degli scores (standardizzati)
plot(Yst(:,1),Yst(:,2),'o')
axislim=axis;
xlabel('Prima componente principale stand. (Indice di povertà)');
ylabel('Seconda componente principale stand. (Indice di malessere delle aziende)');
text(Yst(:,1),Yst(:,2),Xtable.Properties.RowNames)
% 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');

%% Biplot
% Rappresentazione simultanea dei punti riga tramite gli score
% standardizzati e punti colonna tramite le correlazioni tra le variabili
% originarie e le componenti principali
close all
% Rappresentazione nel piano cartesiano dei punti riga tramite gli scores
% standardizzati
plot(Yst(:,1),Yst(:,2),'o')
xlabel('Prima componente principale stand. (Indice di povertà)');
ylabel('Seconda componente principale stand. (Indice di malessere delle aziende)');
text(Yst(:,1),Yst(:,2),Xtable.Properties.RowNames)
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');

% Aggiunta delle frecce che rappresentano le correlazioni tra le variabili
% originarie e le componenti principali
hold('on')
quiver(arx(1,:)',ary(1,:)',arx(2,:)',ary(2,:)')
% viene aggiunto il testo
text(MatrComp(:,1)+dx,MatrComp(:,2)+dy,nameXvars,'Color','Blue');

% axis tight= Fit the axes box tightly around the data by setting the
% axis limits equal to the range of the data.
axis tight