%% LABORATORIO 8 - METODI STATISTICI PER LA BIOINGEGNERIA
%% Canale B - A.A. 2024/2025
%% Docente: Enrico Longato

%% Esercizio 2
clc        
clear all  
close all  

%% Caricamento dati
load bodyfat

%% Definizione dei due problemi di regressione lineare dati
% Cerco gli indici per costruire Y (sempre la stessa) e X, X1
i_y_mask = strcmp('BodyFat', labels);
i_X_mask = ~i_y_mask; % Tutte le altre variabili
i_X1 = [2 3 4]; % Le variabili age, weight e height (indicizzo a mano, per fare prima)

Y = data(:, i_y_mask); 
X_no_intercept = data(:, i_X_mask);
X = [X_no_intercept ones(length(Y), 1)]; % Aggiungo l'intercetta

X1_no_intercept = data(:, i_X1);
X1 = [X1_no_intercept ones(length(Y), 1)]; % Aggiungo l'intercetta

%% Regressioni lineari
% Utilizzo la formula risolutiva per trovare i parametri
beta_hat = (X'*X)\X'*Y;
beta_hat1 = (X1'*X1)\X1'*Y;

% Predizione dei modelli
Y_hat = X*beta_hat;
Y_hat1 = X1*beta_hat1;

% Calcolo dei residui
residuals = Y - Y_hat;
residuals1 = Y - Y_hat1; % NB: è sempre la stessa Y "vera"; cambia la predizione!

% Colleziono le informazioni utili al calcolo di R2, AIC e BIC
SST = sum((Y - mean(Y)).^2); % Sempre la stessa, dipende da Y vero
SSE = sum(residuals.^2);
SSE1 = sum(residuals1.^2);

N = size(X, 1); % Sempre lo stesso
M = size(X, 2); % Numero di parametri intercetta inclusa
M1 = size(X1, 2); % Numero di parametri intercetta inclusa

% Calcolo AIC e BIC
aic = N*log(SSE/N) + 2*M;
aic1 = N*log(SSE1/N) + 2*M1;

bic = N*log(SSE/N) + M*log(N);
bic1 = N*log(SSE1/N) + M1*log(N);

% Calcolo R2 e R2 adjusted
R2 = 1 - SSE/SST;
R21 = 1 - SSE1/SST; % Scusate il nome della variabile; mi sono accorto dopo

% Rispetto alla formula nelle slide di teoria, ho "messo i gradi di libertà
% dove è più facile ricordarsi che vanno", cioè al denominatore della
% quantità da stimare in maniera non polarizzata.
% Nelle slide è stata svolta la "frazione a 4 piani".
R2_adj = 1 - ((SSE/(N-M)) / (SST/(N-1)));
R2_adj1 = 1 - ((SSE1/(N-M1)) / (SST/(N-1)));

% Scrivo e commento i risultati
disp(' ')
disp('Modello con tutte le variabili vs. modello ridotto con le variabili che seguono')
disp(['     - ' strjoin(labels(i_X1), ', ') ' + intercetta'])
disp(['AIC: ' num2str(aic) ' vs. ' num2str(aic1)])
disp(['BIC: ' num2str(bic) ' vs. ' num2str(bic1)])
disp(['R2: ' num2str(R2) ' vs. ' num2str(R21)])
disp(['R2 adjusted: ' num2str(R2_adj) ' vs. ' num2str(R2_adj1)])

disp(' ')
disp(['La riduzione del modello è stata eccessiva: le sole variabili ' ...
    'considerate non sembrano in grado di portare alla definizione di un' ...
    ' modello parsimonioso adeguato.'])