%% 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(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))]; 

% 4. Stima dei parametri beta del modello
beta_hat = 
disp('Beta stimati')
disp(beta_hat)

% 5. Predizione del modello
y_hat =

% 6. Indice di determinazione tra variabile spiegata e predizione del modello
R2 = 
disp('Indice di determinazione R2')
disp(R2)

% 7. Calcolo dei residui
res = 

% 8. Standard error delle misure
sigma_hat_2 =
disp('Stima della varianza dell''errore di misura su hear rate')
disp(sigma_hat_2)

SE =


disp('Standard Error delle misure a disposizione ')
disp(SE)

% 9. Coefficiente di variazione
CV = 
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();
Y_hat_plane = 
surf(X1,X2, Y_hat_plane)

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