clc 
clear all
close all

load data_pca

% I dati sono contenuti nella matrice data. Data contiene 392 osservazioni 
% e 7 variabili, i cui nomi sono in labels e le unità di misura in units.

m = size(data,2);  % numero di variabili
n = size(data,1);  % numero di osservazioni

% Standardizzo le variabili
X = data;
for k =1:m
    X(:,k) = (X(:,k)- mean(X(:,k)))/std(X(:,k));
end
% X è la matrice dei dati standardizzati

% Matrice di covarianza di X
S = cov(X);

% Matrice di correlazione dei dati non standardizzati
R = corr(data);

% Calcolo di autovalori e autovettori di S
[V,D] = eig(S);
% Vettore degli autovalori
lambda = diag(D);
% Ordino gli autovalori in ordine decrescente
[lambda_ord, ind] = sort(lambda,'descend'); 
% Ordino gli autovettori con lo stesso ordine degli autovalori
V_ord = V(:,ind); 

% Applico la trasformazione lineare per passare nello spazio delle
% componenti principali
Y = X*V_ord;

%% Scelta del numero di componenti principali

% Scree plot
figure
plot(lambda_ord,'b-o')
xlabel('# componente principale')
ylabel('autovalore - varianza spiegata')

% Plot della frazione di varianza spiegata dalle prime k componenti
var_tot = sum(lambda); 
var_cum = zeros(1,m); 
for k = 1:m
    var_cum(k) = sum(lambda_ord(1:k));
end
f_var_cum = var_cum/var_tot; 

figure
plot(f_var_cum,'b-o')
xlabel('p componenti principali')
ylabel('Frazione varianza spiegata dalla prime k componenti principali')

% Ricostruzione dei dati originali
X1 = Y*V_ord';
% Verificato che X1 coincide con X

% Sulla base dei plot precedenti scegliamo p = 4 componenti principali.
% Queste consentono di spiegare più del 90% della varianza dei dati
% originali.
p = 4; 

% Dataset trasformato e compresso
Yp = Y(:, 1:p);
% Autovettori che corrispondono alle prime p componenti principali
Vp = V_ord(:,1:p); 
% Ricostruzione dei dati originali usando solo le prime p componenti
% principali
Xr = Yp*Vp';

%% PCA utilizzando la funzione di Matlab

[COEFF, SCORE, LATENT] = pca(X);
% COEFF corrisponde alla matrice V_ord
% SCORE corrisponde alla matrice Y dei dati trasformati
% In LATENT troviamo gli autovalori ordinati dal maggiore al minore
% Nota: per questo dataset pca ci ha restituito una matrice di autovettori
% (COEFF) che differisce dalla nostra V_ord per il segno. Di conseguenza
% anche SCORE differisce da Y per il segno. Il segno non ha importanza
% perché, anche cambiando il segno, gli autovettori rimangono autovettori e
% mantengono le stesse proprietà necessarie a realizzare la PCA. Quindi la
% funzione pca di Matlab ha realizzato una trasformazione lineare dei dati
% equivalente alla nostra che soddisfa le proprietà richieste dalla pca
% tanto quanto la nostra. 

%% PCA usando la SVD

[U D V_svd] = svd(X);
% V_svd a meno del segno coincide con V_ord e COEFF

% Estraggo gli autovalori dalla matrice D
d = diag(D);
d2 = eig(X'*X);
lambda_svd = d.^2/(n-1);
% lambda_svd coincide con lambda_ord e con LATENT