%% ESERCIZIO I 
% Generare 100 osservazioni da una distribuzione normale trivariata con
% vettore delle medie uguale a zero e matrice di covarianze (Sigma)
% riportata di seguito
%          1.00          0.87          0.92
%          0.87          1.00          0.90
%          0.92          0.90          1.00
% Per replicabilità dei risultati utilizzare rng(11).
% Applicare a questa matrice l'analisi in componenti principali partendo
% dalla matrice degli scostamenti dalla media Xtilde. 
% Motivare come mai in questo esempio non si parte dalla matrice degli
% scostamenti standardizzati.
% Mostrare la tabella della varianza spiegata (relativa e cumulate) dalle
% diverse componenti principali. Discutere
% 1) il numero di PC ottimo;
% 2) l'interpretazione della prima PC nel caso in cui le tre variabili
% originarie fossero state collegate a tre test di intelligenza;
% 3) gli angoli tra le frecce nel biplot
%
% Calcolare la matrice di covarianze in maniera matriciale.
% Denominare questa matrice S
% Estrarre dalla matrice S gli autovalori (ordinati in senso non
% decrescente) i rispettivi autovettori e i corrispondenti valori singolari.
% Verificare che il primo valore singolare sia uguale a quello ottenuto
% utilizzando la funzione pcaFS.
% Inserire i valori singolari della matrice S in un vettore colonna di
% lunghezza 3 denominato gamma.
%
% Calcolare la matrice che rappresenta la miglior rappresentazione di rango
% 1 della matrice originaria in termini di scostamenti dalla media.
% Denominare questa matrice Xhat.
% Verificare tramite assert che la matrice Xhat abbia rango 1.
% Qual è la sommma dei quadrati delle differenze tra la matrice Xtilde e la
% matrice Xhat?
% Calcolare la somma dei quadrati degli scostamenti in due modi diversi e
% verificare che entrambi i modi producono gli stessi risultati.
%
% Chiudere tutte le finestre grafiche Costruire il diagramma di dispersione
% 3D tra le 3 variabili originarie in termini di scostamenti dalla media.
% Aggiungere a questo diagramma di dispersione la linea (retta) di massima
% variabilità (linea associata alla prima componente principale).
% Aggiungere al grafico 3D le etichette degli assi ed impostare l'opzione
% di uguale scala tra i diversi assi cartesiani.
%
% Creare una nuova figura con le informazioni della figura precedente e
% l'aggiunta delle proiezioni ortogonali lunga la retta principale
%
% Qual è la lughezza della proiezione ortogonale dell'unità 5 lungo la retta
% principale? Calcolare la distanza euclidea dall'origine dell'unità 5
% proiettata lungo la retta principale.
%
% Calcolare la lunghezza dei semiassi principali dell'ellissoide
% (x-xmedio)' S^-1 (x-xmedio) =1 
%
% Discutere la relazione tra lo scostamento quadratico medio della prima
% componente principali e la lunghezza del primo semiasse principale.

%
% Rappresentare i dati originali tramite coordinate parallele utilizzando
% la funzione parallelcoords. Che tipo di rappresentazione ci si attende?
%
%
% Applicare alla matrice dei dati originaria, il metodo delle k medie
% imponendo 3 gruppi. Utilizzare un numero di repliche pari a 50.
% Visualizzare l'assegnazione delle unità ai diversi gruppi tramite la
% funzione spmplot. Commentare i gruppi che sono stati ottenuti. Scegliere
% automaticamente il numero ottimo dei gruppi tramite la funzione
% evalclusters utilizzando come criterio l'indice di Calinski e Harabasz.
% Utilizzare un numero di gruppi che varia da 1 a 6. Rappresentare
% graficamente l'ouptut della funzione evalclusters.
% Commentare i gruppi che vengono trovati con il valore di k ottimo.
%

% ESERCIZIO II
% Caricare la tabella di contingenza presente dentro il file car.mat di
% FSDA. Questa tabella di contingenza si riferisce alle tipologia di auto
% (righe) e agli attributi che le contraddistinguono (colonne). Calcolare
% il p-value dell'indice Chi2 e discutere la sua relazione con l'inerzia
% totale della tabella di contingenza. Che caratteristiche presentano le
% auto che si trovano nel primo quadrante nel grafico di analisi delle
% corrispondenze?
%
% Utilizzando la funzione CorAnaplot costruire il grafico di analisi delle
% corrispondenze in modo tale che il colore delle etichette dei punti riga
% sia proporzionale alla comunalità (quota di inerzia spiegata) di ogni punto.
% Suggerimento: utilizzare l'opzione  plots.ColorMapLabelRows.
% Commentare il grafico ottenuto. 
% Discutere la differenza in termini di inerzia spiegata tra la Ferrari e
% la Porsche. 
% Calcolare la quota esatta dell'inerzia di Ferrari e Porsche spiegata
% dalle due dimensioni latenti.


%% ESERCIZIO III
% Creare una nuova funzione (denominata erone.m) che dati come argomenti di
% input la lunghezza dei 3 lati di un triangolo, applichi la formula di
% Erone per calcolare l'area del triangolo.
% https://it.wikipedia.org/wiki/Formula_di_Erone. 
% Ad esempio la chiamata a erone(3,4,5) deve produrre come risultato 6.
% Inserire un controllo che fornisca un messaggio di errore se l'utente
% fornisce un numero inferiore a 3 di argomenti di input. Fare in modo che
% la funzione restituisca i due modi per calcolare la formula di Erone
% descritta alla pagina di wikipedia.


%% Generazione della matrice di dimensioni n x 3
rng(11);
Sigma=[1 0.87 0.92;
   0.87 1 0.9;
   0.92 0.9 1];

mu=zeros(3,1);
n=100;
% X è la matrice dei dati che contiene un campione di 100 realizzazione
% dalla distribuzione normale trivariata
X=mvnrnd(mu,Sigma,n);

%% Analisi in componenti principali

% Motivare come mai in questo esempio non si parte dalla matrice degli
% scostamenti standardizzati.
% In questo esempio le variabili presentono tutte lo stesso ordine di
% grandezza (e la stessa unità di misura) di conseguenza non è necessario
% operare la standardizzazione.
Xtable=array2table(X);
outPCA=pcaFS(Xtable,'standardize',false);
% Osservazione: potevo utilizzare anche il comando pcaFS(X).
% Se si utilizza Xtable, le variabili vengono denominate X1, X2 e X3. 

% Mostrare la tabella della varianza spiegata dalle diverse componenti
% principali (PC)
disp('Tabella della varianza spiegata')
disp(outPCA.explainedT)

% Vogliamo prendere un numero di componenti principali che spieghi una
% quota di varianza almeno uguale a 0.95^p=0.95^3=0.8574.
% La prima PC spiega oltre il 93 per cento della variabilità complessiva di
% conseguenza in questo esempio ha senso considerare solo la prima PC.

% Le 3 variabili originarie presentano una correlazione molto forte. La
% prima PC presenta correlazioni positive elevate con le 3 variabili
% originarie.
% Se le tre variabili originarie fossero state collegate a tre test di
% intelligenza, gli scores della prima PC avrebbero permesso di ordinare
% gli individui dal più intelligente al meno intelligente. 
% 
% Commento al biplot: l'angolo tra i vettori X1 e X3 è molto piccolo indica
% che queste due variabili sono correlate tra di loro in maniera forte e
% diretta. La correlazione campionaria tra queste due variabili infatti era
% pari a 0.9329 (la correlazione nell'universo era pari a 0.92).
% L'angolo tra X1 e X3 con X2 è modesto e questo segnala una correlazione
% sempre diretta non modesta tra X1 e X3 con X2. Nella matrice di
% correlazione campionaria r(X1,X2)=0.88 e r(X2,X2)=0.89
%

% meaX = vettore riga delle medie aritmetiche
meaX=mean(X);
% Xtilde matrice degli scostamenti dalla media
Xtilde=X-meaX;
% S = matrice di covarianze calcolata in maniera matriciale.
S=Xtilde'*Xtilde/(n-1);


% Autovalori ed autovettori di S
[Vini,Lambdaini]=eig(S);

[~,ord]=sort(diag(Lambdaini),'descend');
% Lambda contiene sulla diagonale, gli
% autovalori ordinati in senso decrescente
Lambda=Lambdaini(ord,ord);
% V contiene i corrispondenti autovettori
V=Vini(:,ord);
% Verifico che Lambda(1,1) sia uguale a quello che risultava applicando
% direttamente la funzione pcaFS.
assert(abs(Lambda(1,1)-outPCA.explained(1,1))<1e-10,"Errore di programmazione " + ...
    "il primo autovalore è diverso")
% Estrazione dei valori singolari in un vettore colonna denominato gamma
gamma=sqrt(diag(Lambda));


%% Miglior matrice di rango 1
% Calcolare la matrice che rappresenta la miglior rappresentazione di rango
% 1 della matrice originaria in termini di scostamenti dalla media.
% Denominare questa matrice Xhat
% v1 è il primo autovettore (associato al più grande autovalore della
% matrice S).
v1=V(:,1);
Xhat=Xtilde*(v1*v1');
% Verificare tramite assert che la matrice Xhat abbia rango 1.
assert(rank(Xhat)==1,"Attenzione la matrice Xhat non è di rango 1")

% Qual è la sommma dei quadrati delle differenze tra la matrice Xtilde e la
% matrice Xhat?
% sumE= somma dei quadrati degli scostamenti tra la matrice Xtilde e la
% matrice Xhat
sumE=sum((Xtilde-Xhat).^2,"all");
% La somma dei quadrati degli scostamenti tra la matrice Xtilde e la
% matrice Xhat è pari a (n-1) *(secondo autovalore+terzo autovalore)
sumEchk=(n-1)*sum(diag(Lambda(2:end,2:end)));
% Verifico che sumE e sumEchk diano esattamente lo stesso risultato
assert(abs(sumE-sumEchk)<1e-12,"Attenzione la somma dei quadrati " + ...
    "delle differenze è diversa")


% Analisi in componenti principali
% I punti vengono proiettati in maniera ortogonale sulla retta principale.
% Retta principale= retta associata alla direzione di massima variabilità

%% Rappresentare i dati tramite un diagramma di dispersione 3D
% ed aggiungere la retta principale
% Per ulteriori informazioni vedere il file pdf "data science con MATLAB"
close all
scatter3(Xtilde(:,1),Xtilde(:,2),Xtilde(:,3))
hold('on')

% Tutti i punti della matrice Xhat= Xtilde*v1*v1' si trovano lungo una retta
% Prendo due punti qualsiasi per disegnare questa retta.
% Osservazione potevo anche prendere il minimo ed il massimo di Xhat(:,2)
% oppure di Xhat(:,3)
[~,indminX1]=min(Xhat(:,1));
[~,indmaxX1]=max(Xhat(:,1));

% Aggiungo la linea di massima variabilità (linea associata alla prima
% componente principale)
% Affinchè la riga non sia troppo corta faccio riferimento al massimo e
% minimo per una coordinata.

line([Xhat(indminX1,1);Xhat(indmaxX1,1)],[Xhat(indminX1,2);Xhat(indmaxX1,2)],...
    [Xhat(indminX1,3);Xhat(indmaxX1,3)])
% Ad esempio [Xhat(indminX1,1);Xhat(indmaxX1,1)] contiene le coordinate di
% inizio e fine della retta lungo la prima dimensione

% Aggiungere al grafico 3D le etichette degli assi
xlabel('X1')
ylabel('X2')
zlabel('X3')
% Stessa unità di misura per i 3 assi cartesiani
axis equal
% Aggiunta del titolo
title("Diagramma di dispersione 3D con la retta principale")

%% Creare una nuova figura con le informazioni della figura precedente e
% l'aggiunta delle proiezioni ortogonali lunga la retta principale
figure
scatter3(Xtilde(:,1),Xtilde(:,2),Xtilde(:,3))
hold('on')
plot3([Xhat(:,1) Xtilde(:,1)]',[Xhat(:,2) Xtilde(:,2)]',[Xhat(:,3) Xtilde(:,3)]','r')
xlabel('x1')
ylabel('x2')
zlabel('x3')
axis equal
line([Xhat(indminX1,1);Xhat(indmaxX1,1)],[Xhat(indminX1,2);Xhat(indmaxX1,2)],...
    [Xhat(indminX1,3);Xhat(indmaxX1,3)])

% Aggiunta delle proiezioni ortogonali sulla retta principale.
% Aggiungo le linee che partono da Xtilde e arrivano fino a Xhat.
plot3([Xhat(:,1) Xtilde(:,1)]',[Xhat(:,2) Xtilde(:,2)]',[Xhat(:,3) Xtilde(:,3)]','r')
title("Diagramma di dispersione 3D, retta principale e proiezioni " + ...
    "ortogonali lungo la retta principale")
% Provate a fare la rotazione 3D di questo grafico

%% Qual è la lughezza della proiezione ortogonale dell'unità 5 lungo la retta
% principale? Calcolare la distanza euclidea dell'origine dell'unità 5
% rispetto al punto di coordinate 0, 0, 0

% Il vettore |y1|=|Xtilde*v1| contiene le lunghezze delle proiezioni
% ortogonali lungo la retta principale (v. file datascience con MATLAB). 
y1=Xtilde*v1;
lunghezza_proiezione=abs(y1(5));
disp(['La lunghezza della proiezione ortogonale è ' num2str(lunghezza_proiezione)])

% La distanza Euclidea dell'unità 5 proiettata lungo la retta principale
% dall'origine degli assi è esattamente uguale alla proiezione ortogonale
% della riga 5 della matrice Xtilde lungo la retta principale
dist5=sqrt(sum(Xhat(5,:).^2));
disp(['La distanza Euclidea dell''unità 5 dall''origine è ' num2str(dist5)])
% Ovviamente le due quantità coincidono.

% Calcolare la lunghezza dei semiassi principali dell'ellissoide
% (x-xmedio)' S^-1 (x-xmedio) con x=(X1 X2 X3)'
% Le lunghezze dei semiassi principali di questo ellissoide non sono altro
% che i valori singolari della matrice S.
disp("Lunghezza dei semiassi principali dell'ellissoide")
disp(gamma)
% Discutere la relazione tra lo scostamento quadratico medio della prima
% componente principali e la lunghezza del primo semiasse principale. Lo
% scostamento quadratico medio della prima PC è esattamente uguale alla
% lunghezza del primo semiasse principale dell'ellissoide i cui assi
% principali sono ortgonali e passano attraverso le direzioni di massima
% variabilità.

disp(['Scostamento quadratico medio della prima PC ' num2str(std(y1))])
disp(['Lunghezza semiasse principale dell''ellissoide ' num2str(gamma(1))])

%% Rappresentazione in termini di coordinate parallele
figure
% Rappresentare i dati originali tramite coordinate parallele utilizzando
% la funzione parallelcoords. Che tipo di rappresentazione ci si attende.
parallelcoords(X)
% Le variabili sono fortemente correlate in maniera positive di conseguenza
% ci attendiamo poche intersezioni tra le spezzate.

%% Clustering tramite le k-medie
% Applicare il metodo delle k medie imponendo 3 gruppi.
% Utilizzare un numero di repliche pari a 50
% Si inizializza la procedura iterativa facendo
% riferimento a 50 centroidi iniziali.
idx=kmeans(X,3,'Replicates',50);
spmplot(X,idx)
% Commentare i gruppi che sono stati ottenuti. Il metodo delle k-medie
% tende a trovare gruppi di forma sferica. In questo caso, dato che le 3
% variabili sono fortemente correlate, kmeans suddivide il dataset in 3
% gruppi con le seguenti caratteristiche: 
% Un gruppo con i valori elevati di tutte e tre le variabil.i 
% Un gruppo con i valori intermedi di tutte e tre le variabili.
% Un gruppo con i valori bassi di tutte e tre le variabili.

%% Scelta automatica del numero dei gruppi tramite la funzione evalclusters.
% Utilizzare un numero di gruppi che varia da 1 a 6
figure
outcla=evalclusters(X,'kmeans','CalinskiHarabasz','KList',1:6);
plot(outcla)
% Il valore ottimale secondo la procedura automatica, utilizzando una
% sequenza (KList) da 1 a 6 gruppi, è pari a 6.
% idxbest classificazione in presenza del valore di k ottimale.
idxbest=outcla.OptimalY;
spmplot(X,idxbest)
% I gruppi ottenuti presentano forma approssimativamente sferica e valori
% via via crescenti delle tre variabili originarie.

%% Caricare la tabella di contingenza presente dentro il file car.mat di FSDA
clear 
close all
load car
out=CorAna(car);
% Che caratteristiche presentano le auto che si trovano nel primo
% quadrante del grafico di analisi delle corrispondenze?

% Le auto che si trovano nel primo quadrante sono caratterizzate dagli
% attributi FuelEconomy e Value. Questo è il segmento di mercato delle auto
% a bassa cilindrata che consumano poco.

% Calcolare il pvalue dell'indice Chi2 e discutere la sua relazione con
% l'inerzia totale della tabella di contingenza.

% Estraggo il valore del test Chi2 dall'output della funzione CorAna
TestChi2=out.Chi2stat;
% La distribuzione di riferimento del test Chi2 è una v.c. con (I-1)*(J-1)
% gradi di libertà dove I e J sono rispettivamente il numero di righe e di
% colonne della tabella di contingenza.
% Il p-value è la prob. in una v.c. Chi quadrato con (I-1)*(J-1) gradi di
% libertà di ottenere un valore superiore a quello osservato
pval=chi2cdf(TestChi2,(out.I-1)*(out.J-1),'upper');
% Il pvalue è virtualmente uguale a 0. Esiste una relazione significativa
% tra il tipo di auto e gli attributi che la contraddistinguono.
%
% Il valore del test Chi2/n non è altro che l'inerzia totale (varianza
% totale della tabella di contingenza).


% Utilizzando la funzione CorAnaplot costruire il grafico di analisi delle
% corrispondenze in modo tale che il colore delle etichette dei punti riga
% sia proporzionale alla comunalità (quota di inserzia spiegata) di ogni punto.
% Commentare il grafico ottenuto.
plots=struct;
plots.ColorMapLabelRows='CntrbDim2In';
CorAnaplot(out,'plots',plots)
% I punti riga che presentano valori di poco superiori a zero per la prima
% dimensione latente e valori bassi (in modulo) per la seconda dimensione
% latente, tendono ad avere un'inerzia spiegata molto più bassa di quella
% degli altri punti.

% Discutere la differenza in termini di inerzia spiegata tra la Ferrari e
% la Porsche. 
% Ferrari e Porsche si trovano molto vicini nel grafico di analisi delle
% corrispondenzze. Queste due auto sono caratterizzate in larga misura
% dall'attributo style. Mentre la comunalità della Ferrari è alta, quella
% della Porsche è più bassa. In altri termini mentre l'inerzia (distanza
% della Ferrari dal centroide) è spiegata bene dalle prime due dimensioni
% latenti, quella della Porsche è spiegata meno bene.
%
% Calcolare la quota esatta dell'inerzia di Ferrari e Porsche spiegata
% dalle due dimensioni latenti.
out.OverviewRows(["Ferrari","Porsche"],["CntrbDim2In_1","CntrbDim2In_2"])

% Quota di inerzia spiegata dalle prime due dimnsioni latenti di
% "Ferrari" e "Porsche"
disp(sum(out.OverviewRows{["Ferrari","Porsche"],["CntrbDim2In_1","CntrbDim2In_2"]},2))


%% ESERCIZIO III
% Creare una nuova funzione (denominata erone.m) che dati come argomenti di input la lunghezza
% dei 3 lati di un triangolo, tramite la formula di Erone
% https://it.wikipedia.org/wiki/Formula_di_Erone.
% Ad esempio la chiamata a 
% erone(3,4,5) deve produrre come risultato 6.
% Inserire un controllo che fornisca un messaggio di errore se l'utente
% fornisce un numero inferiore a 3 di argomenti di input.
% Fare in modo che la funzione restituisce i due modi per calcolare la
% formula di Erone descritta alla pagina di wikipedia.

% Osservazione il codice che segue deve essere inserito in un file
% denominato erone.m

function [A,Achk]= erone(a,b,c)
%% Formula di erone

if nargin<3
    error('Mi hai dato solo un input, ho bisogno di 3 valori di input')
end

 p=((a+b+c)/2);

 A=sqrt(p*(p-a)*(p-b)*(p-c));

 Achk=sqrt((a+b+c)*(-a+b+c)*(a-b+c)*(a+b-c))/4;

 end