%% OPERAZIONI CON LE MATRICI
% TESTO

% 1a) Creare una matrice identità di dimensione p (assegnare alla variabile p
% un valore qualsiasi) tramite un doppio ciclo for. Creare la matrice
% identità con il comanda eye.
% 1b) Creare la matrice diagonale con elementi sulla diagonale principale
% pari a [2 5 7] e mostrarla nella command
%
% 2) date 3 matrici A, B e C di dimensione mxn, nxo, oxp, verificare la
% regola della moltiplicazione di matrici trasposte   
% (A*B*C)' = C' * B' * A'. Assegnare a m, n, o, p numeri naturali a
% piacere.
%  
%
% 3) Risolvere il sistema di equazioni lineari (A x = b) che segue:
% 3x1+5x2+4x3=25
% 2x1+10x2+x3=25
% 2x1+2x2+9x3=33
%
% 4) dato un vettore colonna x di numeri casuali estratti dalla
% distribuzione normale, scrivere in maniera matriciale la somma dei
% quadrati di x
%
% 5a) dato  un vettore colonna x=(1 2 3)' e un vettore lambda = (2, 3, 6)
% scrivere in maniera matriciale la somma dei quadrati ponderata, ossia
% $\sum_{i=1}^3 x_i^2 \lambda_i$
% 5b) Scrivere la somma dei quadrati ponderati in maniera simbolica


% 6) Generare una matrice di numeri causali dalla distribuzione normale di
% dimensione nxp. Denominare questa matrice X. Assegnare ad n e p numeri
% interi a piacere (con n>p). Calcolare la matrice di covarianze S.
% Verificare empiricamente che la tramite 10000 simulazioni del vettore x
% che la forma quadratica $x'Sx$ è definita positiva. Fare lo storing dei
% 10000 scalari $x'Sx$ in un vettore di dimensione 10000x1 denominato
% formaquad. Mostrare il valore più piccolo di x''Sx delle nsimul
% simulazioni nella command window
%
% Generare una matrice di numeri casuali di dimensione nxp, denominata A,
% dalla distribuzione normale. Premoltiplicare e postmoltiplicare la
% matrice precedente in modo appropriato per poter estrarre l'elemento che
% si trova all'incrocio della riga i e della colonna j (porre i=2, j=5,
% n=10 e p=6)
%
% 8) Data una matrice di numeri casuali di dimensione nxp (ad esempio n=7 e
% p=3), calcolare in maniera matriciale la matrice Xtilde (degli
% scostamenti dalla media) passando attraverso la matrice simmetrica e
% idempotente H
%
% 9) Calcolare in maniera matriciale la matrice di covarianze 
% Confrontare il risultato con quello ottenuto tramite la
% funzione cov
%
% 10) Calcolare in maniera matriciale la matrice degli scostamenti
% standardizzati. Confrontare il risultato con quello ottenuto tramite la
% funzione zscore
% 11) Calcolare:
% 11a) la matrice di correlazione partendo dalla matrice Z
% 11b) la matrice di correlazione tramite la funzione corr 
% 11c) in maniera matriciale tramite le matrici D, X,  H 
%
%
%% OSSERVAZIONE: 
% prima di risolvere la parte di seguito fare l'esercizio nel file
% autov.m
%
% 12) Estrarre gli autovalori ed autovettori della matrice di covarianze e
% verificare la scomposizione spettrale della matrice S
%
% 13) Verificare che la matrice S può essere scritta come
% S= \sum_{i=1}^p \lambda_i v_i v_i'
% Ripetere le stesse operazioni sulla matrice di correlazione R
%  Ricostruzione della matrice di correlazione utilizzando solo i primi due
%  autovalori più grandi (chiamare questa matrice Rhat)
%  Mostarre le differenze tra la matrice R originale e la matrice Rhat
%  tramite grafico heatmap
%
% 14) Scomposizione in valori singolari la matrice Xtilde e verificare la
% scomposizione svd. Verificare che la matrice X tilde può essere scritta come
% Xtilde= \sum_{i=1}^r \gamma_i u_i v_i'
% Verifica del rango della coponente i esima di questa somma ossia il rango della
% matrice  \gamma_i u_i v_i'
% Ricostruire la matrice Xtilde utilizzando solo i primi due addendi della somma
% Rappresentare graficamente tramite heatmap la differenza tra Xtilde e
% Xtildecappello formata dai primi due valori singolari e commentare il grafico
%
%  Verificare che se al posto della matrice Xtilde si sostituisce la
%  matrice Xtildecappello formata dai primi due valori singolari la somma
%  dei quadrati dei residui diviso (n-1) è pari alla somma degli
%  autovalori della matrice di covarianze di S escludendo i due più grandi


%% SOLUZIONE 

%% Matrice identità di dimensione p ottenuta tramite doppio ciclo for
% 1a) Creare una matrice identità di dimensione p (assegnare alla variabile p
% un valore qualsiasi) tramite un doppio ciclo for. Creare la matrice
% identità con il comanda eye.
p=5;
id=zeros(p,p);
for i=1:p
    for j=1:p
        if i==j
            id(i,i)=1;
        end
    end
end
% Il modo veloce per creare la matrice identità è 
I=eye(p);

assert(isequal(I,id),"errore di programmazione ident tramite ciclo for diversa da quella ottenuta con eye")

%% Matrice diagonale
% 1b) Creare la matrice diagonale con elementi sulla diagonale principale
% pari a [2 5 7] e mostrarla nella command
D=diag([2 5 7]);
disp('Esempio di matrice diagonale')
disp(D)

%% 2) verifica  regola moltiplicaione matrice trasposta
m=3;
n=5;
o=7;
p=6;
A=randn(m,n);
B=randn(n,o);
C=rand(o,p);
D=A*B*C;

disp(' Verifica che (A*B*C)''-(C''*B''*A'')  = zeros(p,m)')
disp(D'-(C'*B'*A'))


%% 3) Risoluzione sistemi equazioni lineari
% Risolvere il sistema di equazioni lineari (A x = b) che segue:
% 3x1+5x2+4x3=25
% 2x1+10x2+x3=25
% 2x1+2x2+9x3=33
%
% A =matrice dei coefficienti
% b = vettore contenente i termini noti del sistema
A=[ 3 5 4; 2 10 1; 2 2 9];
b=[25; 25; 33];
x=inv(A)*b;
disp('Soluzione del sistema di equazioni lineari')
disp(x)
% Osservazione: computazionalmente il modo più efficiente per risolvere il
% sistema equazioni lineari è tramite l'operatore \
% x=A\b


%% 4) Somma quadrati in maniera matriciale
n=10; 
x=randn(n,1);
% somma dei quadrati in maniera matriciale
disp('Somma dei quadrati elementi del vettore x')
disp(x'*x)

%% 5a) Somma dei quadrati ponderati
% dato  un vettore x=(1 2 3)' e un vettore lambda = (2, 3 6)
% scrivere in maniera matriciale la somma dei quadrati ponderata, ossia
% $\sum_{i=1}^3 x_i^2 \lambda_i $

x=[1; 2; 3];
lambda=[ 2; 3; 6];
Lambda=diag(lambda);

% somma dei quadrati ponderata
disp('somma quadrati ponderati tramite matrici')
disp(x'*Lambda*x)
disp('somma quadrati ponderati senza matrici')
sum((x.^2).*lambda)


%% 5b) Scrivere la somma dei quadrati ponderati in maniera simbolica
% Definisco i simboli x1 x2 x3 assumendo che siano valori reali 
%{
% Si noti che invece di utilizzare le istruzioni
syms('x1','real')
syms('x2', 'real')
syms('x3','real')
syms('lambda1','lambda2','lambda3')
Lambda=diag([lambda1 lambda2 lambda3]);
x=[x1; x2; x3];
%}
% L'istruzione che segue è un modo molto compatto per definire x e lambda
% come vettori colonna di dimensione 3x1 contenente le variabili a valori
% reali x1, x2, x3, lambda1, lambda2, lambda3
syms('x', 'lambda',[3 1],'real')
Lambda=diag(lambda);

disp('Forma quadratica x''*Lambda*x in maniera simbolica')
disp(x'*Lambda*x)


%% Forme quadratiche definite positive 
% 6) Generare una matrice di numeri causali dalla distribuzione normale di
% dimensione nxp. Denominare questa matrice X. Assegnare ad n e p numeri
% interi a piacere (con n>p). Calcolare la matrice di covarianze S.
% Verificare empiricamente che la tramite 10000 simulazioni del vettore x
% che la forma quadratica $x'Sx$ è definita positiva. Fare lo storing dei
% 10000 scalari $x'Sx$ in un vettore di dimensione 10000x1 denominato formaquad.
% Mostrare il valore più piccolo di x''Sx delle nsimul simulazioni nella
% command window

n=100;
p=3;
X=randn(n,p);
S=cov(X);
% nsimul = numero di simulazioni (ripetizioni dell'esperimento)
nsimul=10000;
% inizializzo il vettore formaquad che conterrà in posizione i il risultato di
% x'Sx per la simulazione i
formaquad=zeros(nsimul,1);
for i=1:nsimul
 x=randn(p,1);
 formaquad(i)=x'*S*x;
 
 assert(formaquad(i)>0,"Errore di programmazione la forma quadratica x'Sx è definita positiva")
end
% Calcolo e mostro il minimo di formaquad
minfq=min(formaquad);
disp('Valore più piccolo di x''Sx in nsimul simulazioni')
disp(minfq)

%% 7) Estrazione elementi in maniera matriciale
% Generare una matrice di numeri casuali di dimensione nxp, denominata A,
% dalla distribuzione normale. Premoltiplicare e postmoltiplicare la
% matrice precedente in modo appropriato per poter estrarre l'elemento che
% si trova all'incrocio della riga i e della colonna j (porre i=2, j=5,
% n=10 e p=6)
n=10;
p=6;
A=randn(n,p);

% Estrazione elemento i,j 
i=2;
j=5;
epre=zeros(n,1);
epre(i)=1;
epost=zeros(p,1);
epost(j)=1;
disp(['Estrazione elmento(' num2str(i) ',' num2str(j) ') di A in maniera matriciale'])
disp(epre'*A*epost)
disp('Verifica')
disp(A(i,j))


%% 8) Matrice scostamenti dalla media in maniera matriciale
n=10;
p=5;
X=randn(n,p);
uno=ones(n,1);
H=eye(n)-uno*uno'/n;
Xtilde=H*X;

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



%% 10) Matrice degli scostamenti standardizzati in maniera matriciale
sigmas=sqrt(diag(S));
D=diag(sigmas);
% invD=inv(D);
invD=D^-1;
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)


%% 11) 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*X' * H * X *invD /(n-1);
disp(Rchk1)


%% OSSERVAZIONE: 
% prima di risolvere la parte di seguito fare l'esercizio nel file
% autov.m

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

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


% 14) Scomposizione in valori singolari la matrice Xtilde e verificare la
% scomposizione svd. Verificare che la matrice X tilde può essere scritta come
% Xtilde= \sum_{i=1}^r \gamma_i u_i v_i'
% Verifica del rango della coponente i esima di questa somma ossia il rango della
% matrice  \gamma_i u_i v_i'
% Ricostruire la matrice Xtilde utilizzando solo i primi due addendi della somma
% Rappresentare graficamente tramite heatmap la differenza tra Xtilde e
% Xtildecappello formata dai primi due valori singolari e commentare il grafico
%
%  Verificare che se al posto della matrice Xtilde si sostituisce la
%  matrice Xtildecappello formata dai primi due valori singolari la somma
%  dei quadrati dei residui diviso (n-1) è pari alla somma degli
%  autovalori della matrice di covarianze di S escludendo i due più grandi

%% 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)
% Buona ricostruzione della matrice di partenza di rango 5 con una matrice
% di rango 2 che utilizza i primi due valori singolari

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

% Somma dei quadrati dei residui tramite l'operatore traccia
lambda_da3apchk=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
trace(Lambda(1:p-2,1:p-2))