%%  Caricamento dati dal fie formulemat.xlsx
Xtable=readtable('formulemat.xlsx','Sheet','dati','Range','A1:D17','ReadRowNames',true);
X=Xtable{:,:};
[n,p]=size(X);

%%  Matrice scostamenti dalla media in maniera matriciale
uno=ones(n,1);
% H matrice idempotente e simmetrica
H=eye(n)-uno*uno'/n;
Xtilde=H*X;

%% 8) Matrice di covarianze in maniera matriciale
S=Xtilde'*Xtilde/(n-1);
disp('Matrice di covarianze tramite la matrice Xtilde')
disp(S)

Schk=X'*H*X/(n-1);
disp('Matrice di covarianze ottenuta come X''H X/(n-1)')
disp(Schk)
disp('Matrice di covarianze ottenuta direttamente tramite la funzione cov')
disp(cov(X))



%% Matrice degli scostamenti standardizzati calcolata in maniera matriciale
sigmas=sqrt(diag(S));
D=diag(sigmas);
invD=inv(D); 
Z=H*X*invD;
disp('Matrice degli scostamenti standardizzati in maniera matriciale')
disp(Z)
disp('Matrice degli scostamenti standardizzati tramite la funzione zscore')
Zchk=zscore(X);
disp(Zchk)


%% Matrice di correlazione in maniera matriciale.
sigmas=sqrt(diag(S));
R=Z'*Z/(n-1);
disp('Matrice di correlazione partendo dalla matrice Z')
disp(R)
disp('Matrice di correlzione tramite la funzione corr')
Rchk=corr(X);
disp(Rchk)
disp('Matrice di correlazione tramite le matrici D, X H ')
Rchk1=invD*S*invD ;
disp(Rchk1)

%% Autovalori ed autovettori
% Estraggo gli autovalori ed autovettori della matrice di covarianze
[V,Lambda]=eig(S);

% Verifico la scomposizione spettrale della matrice S
Schk=V*Lambda*V';

%% Verifica che la matrice S può essere scritta come
% S= \sum_{i=1}^p \lambda_i v_i v_i'
Schk1=zeros(p,p);
for i=1:p
    Schk1=Schk1+Lambda(i,i)*V(:,i)*V(:,i)';
end


%% Estraggo gli autovalori ed autovettori della matrice di correlazione
% Osservazione: questa parte non è presente nel file di Excel
% formulemat(out).xlsx
[VR,LambdaR]=eig(R);

% Verifico la scomposizione spettrale della matrice R
Rchk=VR*LambdaR*VR';

%% Verifica che la matrice R può essere scritta come
% R= \sum_{i=1}^p \lambda_i v_i v_i'
Rchk1=zeros(p,p);
for i=1:p
    Rchk1=Rchk1+LambdaR(i,i)*VR(:,i)*VR(:,i)';
end

%% Ricostruzione della matrice di correlazione utilizzando solo i primi due autovalori più grandi
Rhat=zeros(p,p);
for i=(p-1):p
    Rhat=Rhat+LambdaR(i,i)*VR(:,i)*VR(:,i)';
end
disp(R-Rhat)
heatmap(R-Rhat)

% La parte che segue è contenuta anche nel file di Excel
% formulemat(out).xlsx

%% Scomposizione in valori singolari della matrice Xtilde
[U,Gamma,V]=svd(Xtilde,'econ');
Xtildechk=U*Gamma*V';
disp('Verifica n.1 della scomposizione in valori singolari della matrice Xtilde')
disp(max(Xtilde-Xtildechk,[],'all'))

%% Verifica che la matrice X tilde può essere scritta come
% Xtilde= \sum_{i=1}^r \gamma_i u_i v_i'
Xtildechk1=zeros(n,p);
for i=1:p
    Xtildechk1=Xtildechk1+Gamma(i,i)*U(:,i)*V(:,i)';
end
disp('Verifica n.2 della scomposizione in valori singolari della matrice Xtilde')
disp(max(Xtilde-Xtildechk1,[],'all'))

%% Verifica del rango della matrice  \gamma_i u_i v_i'
i=1; % provare anche con i=2 oppure con i=3
rango=rank(Gamma(i,i)*U(:,i)*V(:,i)');
disp(['Il rango della matrice \gamma_i di dimensioni n x p u_i v_i'' è ' num2str(rango)])

%% Ricostruzione approssimata della matrice Xtilde utilizzando solo i primi due addendi della somma
% Nel ciclo che seuge for i=1:2 si utilizzando solo i primi due valori
% singolari
Xhat=zeros(n,p);
for i=1:2
    Xhat=Xhat+Gamma(i,i)*U(:,i)*V(:,i)';
end

%% Rappresentazione grafica tramite heatmap della differenza tra Xtilde e Xhat
heatmap(Xtilde-Xhat)
% Commento le prime due variabili sono ricostruite molto bene dalla matrice
% che utilizza i primi due valori singolari

%% Calcolo somma dei quadrati dei residui
% Somma dei quadrati dei residui tramite funzione somma
lambda3=sum((Xtilde-Xhat).^2,'all')/(n-1);

% Somma dei quadrati dei residui tramite l'operatore traccia
lambda3chk=trace((Xtilde-Xhat)'*(Xtilde-Xhat))/(n-1);

% Verifica che lambda3 è il quadrato del terzo valore singolare diviso (n-1)
disp('Terzo valore singolare della matrice Xtilde diviso (n-1)')
disp(Gamma(3,3)^2/(n-1))

disp('Terzo autovalore della matrice S (autovalore più piccolo)')
% Osservazione nella matrice Lambda gli autovalori sono ordinati dal più
% piccolo al più grande
disp(Lambda(1,1))