clc
clear all
close all

load dati_pressione

x = pressione_sistolica; 
y = pressione_diastolica;
n = length(x);

%% Analisi di correlazione

% Coefficiente di correlazione calcolato usando lo stimatore campionario
r = sum((x-mean(x)).*(y-mean(y)))/((n-1)*std(x)*std(y));
disp(['Coefficiente di correlazione: ' num2str(r)])

figure
scatter(x,y)
xlabel('Pressione sistolica [mmHg]')
ylabel('Pressione diastolica [mmHg]')
box on
set(gca,'fontsize',14)
text(100,110,['r=' num2str(r)],'fontsize',14)

% Coefficiente di correlazione calcolato dalla funzione Matlab
r_matlab = corr(x,y);
disp(['Coefficiente di correlazione calcolato con la funzione corr: ' num2str(r_matlab)])

%% Analisi di regressione semplice
X = [ones(n,1) x];

% Stima dei coefficienti con il metodo dei minimi quadrati lineari
beta = inv(X'*X)*X'*y;

% Predizione del modello
y_hat = X*beta; 

% Residui
residui = y-y_hat; 

% Calcolo di R2
SSE = sum(residui.^2);
SST= sum((y-mean(y)).^2);
R2 = 1-SSE/SST;
disp(['Coefficiente di determinazione: ' num2str(R2)])

%% Confronto dei risultati

% Quadrato del coefficiente di correlazione
disp(['Coefficiente di correlazione al quadrato: ' num2str(r^2)])

% Calcolo il coefficiente di correlazione utilizzando il coeffciente beta della regressione lineare
r_reg = beta(2)*std(X(:,2))/std(y);
disp(['Coefficiente di correlazione calcolato a partire dalla regressione: ' num2str(r_reg)])

%% Analisi di regressione lineare multipla

x1 = pressione_sistolica;
x2 = eta;
X = [ones(n,1) x1 x2];

% Stima dei coefficienti con il metodo dei minimi quadrati lineari
beta = inv(X'*X)*X'*y;

% Predizione del modello
y_hat = X*beta; 

% Residui
residui = y-y_hat; 

% Calcolo di R2
SSE = sum(residui.^2);
SST= sum((y-mean(y)).^2);
R2 = 1-SSE/SST;
disp(' ')
disp('Modello di regressione lineare multipla: ')
disp(['Coefficiente di determinazione: ' num2str(R2)])

% Calcolo del coefficiente di correlazione lineare tra y e y_hat
r = corr(y,y_hat);
disp(['Coefficiente di correlazione tra y e y_hat al quadrato: ' num2str(r^2)])

figure
scatter(y,y_hat)
xlabel('Pressione diastolica reale [mmHg]')
ylabel('Pressione diastolica predetta [mmHg]')
box on
set(gca,'fontsize',14)
text(50,100,['r=' num2str(r)],'fontsize',14)