%% 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 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. 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. 9) Fare il clear della memoria 10) Reimportare in formato table (utilizzare come nomi delle righe le province) i dati presenti nel file benessere.xlsx Effettuare la scomposizione in valori singolari (svd) della matrice degli scostamenti standardizzati. Z= U \Gamma V' 11) Utilizzando le informazioni presenti nella svd. Costruire un biplot in cui: 11a) i punti riga sono rappresentati dalle componenti principali standardizzate e i punti colonna tramite frecce la cui lunghezza rappresenta le correlazioni tra le variabili originarie e le componenti principali; 11b) i punti riga sono rappresentati dalle componenti principali non standardizzate e i punti colonna tramite frecce la cui lunghezza rappresenta gli autovettori; 11c) i punti riga sono rappresentati dalle coordinate della matrice nx2  \sqrt(n-1)^\omega *U_(2)*\Gamma_(2)^\alpha e i punti colonna tramite frecce le cui coordinate sono date dalla matrice \sqrt(n-1)^(1-\omega) *Gamma_(2)^(1-\alpha) *V_(2)^T. con \omega scalare che può assumere valori, 0, 0.1, 0.2, ..., 1 e con \alpha scalare che può assumere valori, 0, 0.1, 0.2, ..., 1 Mostrare che il caseo 11a) si ottiene ponendo \alpha=0 e \omega=1 e che il caso 11b) si ottiene ponendo \alpha=1 e \omega=1 12) Trovare gli autovettori, gli score non standardizzati e gli autovalori della matrice di correlazione utilizzando direttamente la funzione MATLAB pca
%% Soluzione
%% 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);

% varalbs = cell che contiene i nomi delle varaibili
varlabs=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);

% 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);
[~,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
trace(La);
la=diag(La);

%
%% Calcolare la matrice che contiene le
% componenti principali (Y=Z V) e la matrice che contiene le componenti

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

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

%
%% 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];
namerows=cellstr([repmat('PC',p,1) num2str((1:p)')]);
namecols={'Autovalori' 'Var_spiegata' 'Var_cum_spiegata'};
autovaltable=array2table(autoval,'RowNames',namerows,'VariableNames',namecols);

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

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

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

%
%% 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',varlabs,...
    'VariableNames',{'Communalities'});


% 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 bixplot e poi tramite l'implementazione manuale
%% Rappresentazione grafica delle correlazioni tramite la funzione biplot

biplot(MatrComp(:,1:2),'varlabels',varlabs)
%   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 tramite implementazione manuale delle correlazioni

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,varlabs);
% 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 frecce 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,varlabs);

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

% 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 prliminarmente 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,varlabs,'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

%
%% Punto 10)
% Reimportare in formato table (utilizzare come nomi delle righe le
% province) i dati presenti nel file benessere.xlsx
% Effettuare la scomposizione in valori singolari (svd) della matrice degli
% scostamenti standardizzati. Z= U \Gamma V'
% 11) Utilizzando le informazioni presenti nella svd.
% Costruire un biplot in cui:
% 11a) i punti riga sono rappresentati dalle componenti principali
% standardizzate e i punti colonna tramite frecce la cui lunghezza
% rappresenta la correlazioni tra le variabili originarie e le componenti
% principali;
% 11b) i punti riga sono rappresentati dalle componenti principali non
% standardizzate e i punti colonna tramite frecce la cui lunghezza
% rappresenta gli autovettori;
% 11c) i punti riga sono rappresentati dalle coordinate della matrice nx2
% $\sqrt(n-1)^\omega *U_(2)*\Gamma_(2)^\alpha$
% e i punti colonna tramite frecce le cui coordinate sono date dalla
% matrice $\sqrt(n-1)^(1-\omega) *Gamma_(2)^(1-\alpha) *V_(2)^T$.
% con $\omega$ scalare che può assumere valori, 0, 0.1, 0.2, ..., 1
% e con $\alpha$ scalare che può assumere valori, 0, 0.1, 0.2, ..., 1
% Mostrare che il caseo 11a) si ottiene ponendo \alpha=0 e \omega=1 e che
% il caso 11b) si ottiene ponendo \alpha=1 e \omega=1

Xtable=readtable('benessere.xlsx','ReadRowNames',true);
% X = matrice di double senza nomi delle righe e nomi delle colonne
X=table2array(Xtable);
[n,p]=size(X);

% Standardizzo i dati
Z=zscore(X);

[U,Gammastar,V]=svd(Z,'econ');
% Controllo sulla svd
max(max(abs(Z-U*Gammastar*V')))
sqn1=sqrt(n-1);
Gamma=Gammastar/sqn1;

%
%% 11a) soluzione
% i punti riga sono rappresentati dalle componenti principali standardizzate e i punti colonna tramite frecce la cui lunghezza rappresenta la correlazioni tra le variabili originarie e le componenti principali;

close all
hold('on')
plot(sqn1*U(:,1),sqn1*U(:,2),'o')
text(sqn1*U(:,1),sqn1*U(:,2),Xtable.Properties.RowNames)

Vgam=V*Gamma;
zeroes = zeros(p,1); 
quiver(zeroes,zeroes,Vgam(:,1),Vgam(:,2))
varlabs=Xtable.Properties.VariableNames;
dx=0.02;
dy=0.03;
text(Vgam(:,1)+dx,Vgam(:,2)+dy,varlabs,'Color','b');

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');

%
%% 11b) soluzione
% i punti riga sono rappresentati dalle componenti principali non standardizzate e i punti colonna tramite frecce la cui lunghezza rappresenta gli autovettori;

close all
hold('on')
plot(sqn1*U(:,1)*Gamma(1,1),sqn1*U(:,2)*Gamma(2,2),'o')
text(sqn1*U(:,1)*Gamma(1,1),sqn1*U(:,2)*Gamma(2,2),Xtable.Properties.RowNames)

Vgam=V;
zeroes = zeros(p,1); 
quiver(zeroes,zeroes,Vgam(:,1),Vgam(:,2))
varlabs=Xtable.Properties.VariableNames;
dx=0.02;
dy=0.03;
text(Vgam(:,1)+dx,Vgam(:,2)+dy,varlabs,'Color','b');

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');

%
%% 11c) soluzione
% i punti riga sono rappresentati dalle coordinate della matrice nx2  \sqrt(n-1)^\omega *U_(2)*\Gamma_(2)^\alpha e i punti colonna tramite frecce le cui coordinate sono date dalla matrice \sqrt(n-1)^(1-\omega) *Gamma_(2)^(1-\alpha) *V_(2)^T. con \omega scalare che può assumere valori, 0, 0.1, 0.2, ..., 1 e con \alpha scalare che può assumere valori, 0, 0.1, 0.2, ..., 1 Mostrare che il caso 11a) si ottiene ponendo \alpha=0 e \omega=1 e che il caso 11b) si ottiene ponendo \alpha=1 e \omega=1 i punti riga sono rappresentati dalle componenti principali non standardizzate e i punti colonna tramite frecce la cui lunghezza rappresenta gli autovettori;

close all
hold('on')
omega = 0.7;
alpha= 0.3;

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

plot(PuntiRiga(:,1),PuntiRiga(:,2),'o')
text(PuntiRiga(:,1),PuntiRiga(:,2),Xtable.Properties.RowNames)

zeroes = zeros(p,1); 
quiver(zeroes,zeroes,PuntiColonna(:,1),PuntiColonna(:,2))
varlabs=Xtable.Properties.VariableNames;
dx=0.02;
dy=0.03;
text(PuntiColonna(:,1)+dx,PuntiColonna(:,2)+dy,varlabs,'Color','b');

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');

%
%% 12) Trovare gli autovettori, gli score non standardizzati e gli
% autovalori della matrice di correlazione utilizzando direttamente la funzione MATLAB pca

% Analisi in componenti principali tramite la funzione pca
% V=  è la matrice degli autovettori
% score = componenti principali non standardizzate
% latent = autovalori ordinati dal più grande al più piccolo
[V,score,latent]=pca(Z);