%% 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 e 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 prime due
% componenti principali (matrici di componenti).
% Verificare che la somma
% dei quadrati delle colonne della matrice di componenti sono uguali ai
% rispettivi autovalori.
% 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 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 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','Sheet','X (database originale)','ReadRowNames',true,'Range','A3:H106');
% X = matrice di double senza nomi delle righe e nomi delle colonne
X=table2array(Xtable);

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

[n,p]=size(X);

%% 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.

% 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 è una matrice di dimensione px3 che contiene
% Prima colonna: autovalori ordinati
% Seconda colonna: 100*autovalori/p
% Terza colonna: somma cumulata in percentuale
autoval=[lasor 100*(lasor)/p 100*cumsum(lasor)/p];
namerows=cellstr([repmat('PC',p,1) num2str((1:p)')]);
namecols={'Autovalori' 'Var_spiegata' 'Var_cum_spiegata'};
autovaltable=array2table(autoval,'RowNames',namerows,'VariableNames',namecols);


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


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

% 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
% namePCS contiene la sequenza di etichette PC1, ..., PCp
% namePCs è un cell array
% L'istruzione repmat('PC',p,1) ripete per p righe e una colonna la stringa
% PC. In questo caso crea il vettore colonne con il character 'PC' ripetuta
% p volte. 
% num2str((1:p)') crea il vettore colonna di characters
% '1' '2' ..., 'p'
namePCs=cellstr([repmat('PC',p,1) num2str((1:p)')]);

% In alternativa, la sequenza di etichette PC1, ..PCp poteva essere creata
% tramite il concatenamento di stringhe come segue.
namePCS="PC"+(1:p)';
% In questo caso namePCs è uno string array

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.

% Le label sull'asse x del grafico a barre sono i nomi delle variabili, il
% secondo argomento di categorical specifica che vogliano usare l'ordine
% delle labels specificato dentro varlabs e non l'ordine alfabetico
xlabels=categorical(nameXvars,nameXvars);

figure
for j=1:2
    subplot(2,1,j)
    b=bar(xlabels, MatrComp(:,j),'g');
    title(['Correlazioni con PC' num2str(j)])
    % xtips e ytips sono le coordinate numeriche dove inserire le etichette
    % del valore numerico sopra ogni barra
    xtips = b.XEndPoints;
    ytips = b.YEndPoints;
    barlabels = string(round(MatrComp(:,j),2));
    text(xtips,ytips,barlabels,'HorizontalAlignment','center',...
        'VerticalAlignment','bottom')
end



%% 5) Verificare tramite la funzione corr le correlazioni tra la prima
% componente principale 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 dalle 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',{'Comunalità'});
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 delle correlazioni utilizzando la funzione quiver
% Rappresentazione delle frecce si possono tracciare utilizzando la funzione
% quiver
close all
zeroes = zeros(p,1);
% Frecce che partono dall'origine e arrivano fino a MatrComp
quiver(zeroes,zeroes,MatrComp(:,1),MatrComp(:,2))
% Label delle frecce
text(MatrComp(:,1),MatrComp(:,2),nameXvars,...
    'VerticalAlignment','bottom','HorizontalAlignment','center');
% Vengono aggiunti gli assi cartesiani
xline(0)
yline(0)
% Aggiunta delle label sugli assi
xlabel('Prima PC: Indice di povertà');
ylabel('Seconda PC: Indice di malessere delle aziende');
% 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')
xlabel('Prima PC: Indice di povertà');
ylabel('Seconda PC: Indice di malessere delle aziende');
text(Y(:,1),Y(:,2),Xtable.Properties.RowNames)
% Vengono aggiunti gli assi cartesiani
xline(0)
yline(0)

%% Calcolo delle componenti principali standardizzate
% Yst= componenti principali standardizzate cov(Yst) = matrice identità
Yst=Y*sqrt(inv(La));
% Rappresentazione nel piano cartesiano degli scores (standardizzati)
plot(Yst(:,1),Yst(:,2),'o')
xlabel('Prima PC: Indice di povertà');
ylabel('Seconda PC: Indice di malessere delle aziende');
text(Yst(:,1),Yst(:,2),Xtable.Properties.RowNames)
% Vengono aggiunti gli assi cartesiani
xline(0)
yline(0)


%% 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)))))