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

clc
clear all
close all

%% Definizione del seed del generatore di numeri casuali
% Serve perché le parti "casuali" degli algoritmi implementati da MATLAB
% diano sempre lo stesso risultato
rng default
rng(42) % 42 è il seed; su tutti i computer, usando 42, dovrebbe venire lo stesso identico risultato

%% Caricamento dati
load bodyfat_extended

%% Definizione del problema di regressione lineare
% Cerco gli indici per costruire Y e X
i_y_mask = strcmp('BodyFat', labels);
i_X_mask = ~i_y_mask;

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

%% Stima dei beta con fitlm (già visto)
mdl = fitlm(X_no_intercept, Y, 'Intercept', true);
beta_hat = mdl.Coefficients.Estimate;

%% Regolarizzazione L1 (svolto)
% NB: la function lasso implementa, in realtà, il metodo elastic net.
% Quindi, per alpha circa uguale a zero (ma non = 0 perché dà errore), è la
% ridge (L2); per alpha = 1 (questa volta uguale) è la LASSO (L1)
lambda_range = logspace(-3, 3, 50); % "50 numeri logaritmicamente equispaziati da 10^-3 a 10^3"
[beta_hat_L1, info_L1] = lasso(X_no_intercept, Y, 'Intercept', true, ...
    'Standardize', true, ... % Standardizzare prima di regolarizzare è sempre bene
    'Alpha', 1, ... % Forza la lasso
    'lambda', lambda_range, ... % Intervallo GEOMETRICO dei lambda
    'CV', 5); % Numero di iterazioni della cross-validazione

% Dentro a beta_hat_L1 ci sono i valori dei parametri corrispondenti a un
% determinato valore di lambda in lambda_range, uno per colonna
% Le intercette sono in info_L1.Intercept, sempre una per ogni lambda
% Dobbiamo trovare il valore "ottimo", cioè quello che minimizza l'MSE, di
% lambda, ed estrarre i parametri beta_hat corrispondenti
[mse_min_L1, i_lambda_best_L1] = min(info_L1.MSE);
lambda_best_L1 = lambda_range(i_lambda_best_L1);
beta_hat_best_L1 = beta_hat_L1(:, i_lambda_best_L1);
intercept_best_L1 = info_L1.Intercept(i_lambda_best_L1);

%% Regolarizzazione L2 (proposto)
[beta_hat_L2, info_L2] = lasso(X_no_intercept, Y, 'Intercept', true, ...
    'Standardize', true, ... % Standardizzare prima di regolarizzare è sempre bene
    'Alpha', 1e-32, ... % Numero piccolissimo; forza la ridge
    'lambda', lambda_range, ... % Intervallo GEOMETRICO dei lambda
    'CV', 5); % Numero di iterazioni della cross-validazione

% Stesso procedimento di prima
[mse_min_L2, i_lambda_best_L2] = min(info_L2.MSE);
lambda_best_L2 = lambda_range(i_lambda_best_L2);
beta_hat_best_L2 = beta_hat_L2(:, i_lambda_best_L2);
intercept_best_L2 = info_L2.Intercept(i_lambda_best_L2);

%% Confronto tra i tre modelli (non regolarizzato, regolarizzati L1 e L2)
% NB: si poteva fare tutto con un bel ciclo for, ma così è più chiaro

% Recupero tutte le informazioni "necessarie" per le formule qui sotto
N = length(Y); % Sempre lo stesso

M = length(beta_hat); % C'era già l'intercetta
M_L1 = sum(abs(beta_hat_best_L1) > 0) + 1; % Numero di parametri diversi da 0 più l'intercetta che era a parte
M_L2 = M; % Nessun coefficiente può andare a zero

SST = sum((Y - mean(Y)).^2); % Sempre la stessa; dipende dal valore vero di Y

% Attenzione a dove sono i coefficienti beta e l'intercetta in ciascun
% metodo!
residuals = Y - (X_no_intercept*beta_hat(2:end) + beta_hat(1)); 
residuals_L1 = Y - (X_no_intercept*beta_hat_best_L1 + intercept_best_L1);
residuals_L2 = Y - (X_no_intercept*beta_hat_best_L2 + intercept_best_L2);

SSE = sum(residuals.^2);
SSE_L1 = sum(residuals_L1.^2);
SSE_L2 = sum(residuals_L2.^2);

% AIC
aic = N*log(SSE/N) + 2*M;
aic_L1 = N*log(SSE_L1/N) + 2*M_L1;
aic_L2 = N*log(SSE_L2/N) + 2*M_L2;

% BIC
bic = N*log(SSE/N) + log(N)*M;
bic_L1 = N*log(SSE_L1/N) + log(N)*M_L1;
bic_L2 = N*log(SSE_L2/N) + log(N)*M_L2;

% R2_adjusted
R2_adj = 1 - ((SSE/(N-M)) / (SST/(N-1)));
R2_adj_L1 = 1 - ((SSE_L1/(N-M_L1)) / (SST/(N-1)));
R2_adj_L2 = 1 - ((SSE_L2/(N-M_L2)) / (SST/(N-1)));

% Commento
disp(' ')
disp(['AIC: ' num2str(aic) ' vs. ' ...
    num2str(aic_L1) ' (L1) vs. ' ...
    num2str(aic_L2) ' (L2)'])
disp(['BIC: ' num2str(bic) ' vs. ' ...
    num2str(bic_L1) ' (L1) vs. ' ...
    num2str(bic_L2) ' (L2)'])
disp(['R2 adjusted: ' num2str(R2_adj) ' vs. ' ...
    num2str(R2_adj_L1) ' (L1) vs. ' ...
    num2str(R2_adj_L2) ' (L2)'])

disp(' ')
disp('La regolarizzazione L1 porta ad un modello che è, rispetto a quello non regolarizzato:')
disp('  - leggermente migliore secondo AIC')
disp('  - migliore secondo BIC')
disp('  - comparabile secondo R2 adjusted (praticamente identico)')
disp(' ')
disp(['Conclusione: preferiamo il modello regolarizzato L1,' ...
    ' che ha ' num2str(M_L1) ' variabili invece di ' num2str(M) ' a circa parità di performance'])
disp(' ')
disp('La regolarizzazione L2 porta ad un modello che è, rispetto a quello non regolarizzato:')
disp('  - peggiore secondo AIC')
disp('  - peggiore secondo BIC')
disp('  - peggiore secondo R2 adjusted')
disp(' ')
disp(['Conclusione: come ci si poteva aspettare, siccome il numero di variabili non cambia rispetto al modello' ...
    ' non regolarizzato, ma, per costruzione, la regolarizzazione introduce un bias nelle stime' ...
    ' (che ha una sua utilità per scopi di predizione, ma questo fatto non viene tenuto in considerazione da BIC/AIC/R2adjusted), ' ...
    'BIC/AIC/R2adjusted non rilevano alcun vantaggio nell''utilizzare il modello regolarizzato L2'])
disp('NB: la differenza tra "comparabile" e "leggermente peggiore" è una questione di sensibilità personale; all''esame basta giustificare')