%% LABORATORIO 5 - 6 Novembre 2024 - METODI STATISTICI PER LA BIOINGEGNERIA

close all;
clearvars;
clc;

%% ESERCIZIO 1:

%% Parte A - Visualizzazione dei dati

% 1. Caricamento dei dati contenuti in data_ELSA.mat
load('data_Lab5.mat')

% 2. Eliminare dalla matrice elsa eventuali valori mancanti, negativi o non fisiologici
% any(:,2) esegue l'operazione per righe (cioè i soggetti)
flag_delete = any(isnan(elsa) | isinf(elsa) | (elsa<0),2);
elsa(flag_delete,:) = [];

% 3. Estrarre le variabili di interesse: Frequenza cardiaca (hr) (variabile dipendente) e
% Pressione sistolica (bp)
hr = elsa(:,4);
bp = elsa(:,15);

% 4. Realizzare il plot di correlazione (corrplot) delle variabili di interesse
figure(1)
corrplot(table(hr,bp))

%% Parte B - Impostazione del modello e stima dei parametri

% 5. Selezionare la heart rate come variabile dipendente da stimare tramite
% la blood pressure mediante relazione lineare (Modello: Y=X*Beta+q)
Y=(hr);
X=[bp ones(size(bp))];

% 6. Stima dei parametri beta del modello conformula chiusa
beta_hat = (X'*X)\(X')*Y;
disp('Beta stimati')
disp(beta_hat)

% 7. Predizione del modello
y_hat = X*beta_hat;


%% Parte C - Verificare i risultati ottenuti in termini di errori di misura e precisione delle stime 

% 8. Indice di determinazione tra variabile spiegata e predizione del modello
R2 = corr(Y,y_hat)^2;
disp('Indice di determinazione R2')
disp(R2)

% 9. Calcolo dei residui
res = Y-y_hat;

% 10. Varianza dell'errore  
n = length(Y);
m = length(beta_hat);
sigma_hat_2 = (res'*res)/(n-m);

disp('Stima della varianza dell''errore di misura')
disp(sigma_hat_2)


% 11. Standard error delle stime e Coefficiente di variazione
SE =sqrt((sigma_hat_2)*diag(inv(X'*X))); 
disp('Standard Error delle misure')
disp(SE)

CV = SE./abs(beta_hat)*100;
disp('Coefficienti di variazione percentuali dei beta')
disp(CV)

%% Parte D - Visualizzazione grafica dei risultati

% 12. Visualizzare mediante scatter plot Y e y_hat  
figure(3) 
scatter(Y, y_hat,'*')
grid on 
xlabel('Y') 
ylabel('Yhat') 
title(['Scatterplot ' 'Y' '-' 'Yhat'])

% 13. Preparazione del modello utilizzato beta*X + q
x = [linspace(0,220,1000)'  ones(1000,1)];
y_hat2 = x*beta_hat;

% 14. Visualizzare nella stessa figura (subplot): dati e predizione; residui 
figure(4)
subplot(211)
plot(X(:,1), Y,'+')
hold on
plot(x(:,1), y_hat2,'-')
ylim([-30 110])
xlim([0 220])
legend({'X vs Y','Linear model'},'Location','northwest')
xlabel('Pressione sistolica [mmHg]')
ylabel('Frequenza cardiaca [bpm]')

subplot(212)
plot(res)
title('Residui (Y - Yhat)')

%% ESERCIZIO 2

close all;

% 1. Estrarre le variabili di interesse: heart rate (heart_rate) (variabile dipendente),
% blood pressure (systolic_bp) e glucosio (gluc). Concatenare in una
% matrice di due colonne la blood pressure e il glucosio

hr   = elsa(:,4);
bp   = elsa(:,15);
gluc = elsa(:,9);

% 2. Realizzare il plot di correlazione (corrplot) delle variabili di
% interesse
%figure(1)
%corrplot(table(hr,bp,gluc))

% 3. Selezionare la heart rate come variabile dipendente da stimare tramite
% la blood pressure e il glucosio mediante relazione lineare (Modello: Y=X*Beta+q)
Y=hr;
X=[bp gluc ones(size(bp,1),1)];

% 4. Stima dei parametri beta del modello
beta_hat = (X'*X)\(X')*Y;
disp('Beta stimati')
disp(beta_hat)

% 5. Predizione del modello
y_hat = X*beta_hat;

% 6. Indice di determinazione tra variabile spiegata e predizione del modello
R2 = corr(Y,y_hat)^2;
disp('Indice di determinazione R2')
disp(R2)

% 7. Calcolo dei residui
res = Y-y_hat;

% 8. Varianza dell'errore  
sigma_hat_2 = (res'*res)/(length(Y)-length(beta_hat));
disp('Stima della varianza dell''errore di misura su hear rate')
disp(sigma_hat_2)


% 9. Standard error delle stime e Coefficiente di variazione
SE =sqrt((sigma_hat_2)*diag(inv(X'*X)));
disp('Standard Error delle misure a disposizione ')
disp(SE)

CV = SE./abs(beta_hat)*100;
disp('Coefficienti di variazione')
disp(CV)

% 10. Visualizzare mediante scatter plot Y e y_hat  
figure(3) 
scatter(Y, y_hat,'*')
grid on 
xlabel('Y') 
ylabel('Yhat') 
title(['Scatterplot ' 'Y - Yhat'])

% 11. Visualizzare nella stessa figura (subplot): dati e predizione; residui 
figure(4)
subplot(211)
plot3(X(:,1),X(:,2), Y,'+')
xlabel('Pressione sistolica [mmHg]')
zlabel('Frequenza cardiaca [bpm]')
ylabel('Concentrazione di glucosio [mg/dl]')

hold on
% BONUS
[X1,X2] = meshgrid(linspace(80,220,10),linspace(60,290,10));
Y_hat_plane = X1*beta_hat(1)+X2*beta_hat(2)+beta_hat(3);
surf(X1,X2, Y_hat_plane)

subplot(212)
plot(res)
title('Residui (Y - Yhat)')

%% ESERCIZIO 1 BONUS (Z-SCORE):

%% Parte A - Visualizzazione dei dati

% 1. Caricamento dei dati contenuti in data_ELSA.mat
load('data_Lab5.mat')

% 2. Eliminare dalla matrice elsa eventuali valori mancanti, negativi o non fisiologici
flag_delete = any(isnan(elsa) | isinf(elsa) | (elsa<0),2);
elsa(flag_delete,:) = [];

% 2a. Z-SCORE
elsa=zscore(elsa);

% 3. Estrarre le variabili di interesse: Frequenza cardiaca (hr) (variabile dipendente) e
% Pressione sistolica (bp)
hr = elsa(:,4);
bp = elsa(:,15);

% 4. Realizzare il plot di correlazione (corrplot) delle variabili di interesse
figure(1)
corrplot(table(hr,bp))

%% Parte B - Impostazione del modello e stima dei parametri

% 5. Selezionare la heart rate come variabile dipendente da stimare tramite
% la blood pressure mediante relazione lineare (Modello: Y=X*Beta+q)
Y=hr;
X=bp;

% 6. Stima dei parametri beta del modello conformula chiusa
beta_hat = (X'*X)\(X')*Y;
disp('Beta stimati')
disp(beta_hat)

% 7. Predizione del modello
y_hat = X*beta_hat;


%% Parte C - Verificare i risultati ottenuti in termini di errori di misura e precisione delle stime 

% 8. Indice di determinazione tra variabile spiegata e predizione del modello
R2 = corr(Y,y_hat)^2;
disp('Indice di determinazione R2')
disp(R2)

% 9. Calcolo dei residui
res = Y-y_hat;

% 10. Varianza dell'errore
n = length(Y);
m = length(beta_hat);
sigma_hat_2 = (res'*res)/(n-m);

disp('Stima della varianza dell''errore di misura')
disp(sigma_hat_2)

% 11. Standard error delle stime e Coefficiente di variazione
SE =sqrt((sigma_hat_2)*diag(inv(X'*X))); 
disp('Standard Error delle stime dei beta')
disp(SE)

CV = SE./abs(beta_hat)*100;
disp('Coefficienti di variazione percentuali dei beta')
disp(CV)

%% Parte D - Visualizzazione grafica dei risultati

% 12. Visualizzare mediante scatter plot Y e y_hat  
figure(3) 
scatter(Y, y_hat,'*')
grid on 
xlabel('Y') 
ylabel('Yhat') 
title(['Scatterplot ' 'Y' '-' 'Yhat'])

% 13. Preparazione del modello utilizzato beta*X + q
x = linspace(-3,3,1000)';
y_hat2 = x*beta_hat;

% 14. Visualizzare nella stessa figura (subplot): dati e predizione; residui 
figure(4)
subplot(211)
plot(X(:,1), Y,'+')
hold on
plot(x(:,1), y_hat2,'-')
ylim([-3 3])
xlim([-3 3])
legend({'X vs Y','Linear model'},'Location','northwest')
xlabel('Pressione sistolica [mmHg]')
ylabel('Frequenza cardiaca [bpm]')

subplot(212)
plot(res)
title('Residui (Y - Yhat)')