%% LABORATORIO 9 - METODI STATISTICI PER LA BIOINGEGNERIA 
close all;
clearvars;
clc;

%% ESERCIZIO 1: REGRESSIONE E METODI DI SHRINKAGE
% LOADING E PRE-PROCESSING %

% 1. Fissare il seed
rng("default");
rng(1);

% 2. Caricare i dati e chiamare n il numero di soggetti
load('data_lab9.mat');
n=size(data, 1);

% 3. Creare la variabile X contenente le variabili indipendenti (fino a end-1) 
% e Y contenente la variabile dipendente (ultima colonna). 
X = data(:,1:end-1);
Y = data(:,end);


% 4. Fittare i dati con il metodo dei minimi quadrati lineari il modello 
% di regressione, calcolare il coefficiente di determinazione e la conseguente 
% predizione 
B_reg = (X'*X)\X'*Y;
Y_hat_reg=X*B_reg;
R2 = corr(Y, Y_hat_reg).^2;

% 5. Rappresentare in un plot: 1) l'andamento di Y attraverso i soggetti, con sovrapposta la stima del modello 2) scatter plot del dato vs la predizione.
figure(1)
subplot(121)
plot(1:n, Y, 'o--'), hold on, plot(1:n, Y_hat_reg, 'r', 'linewidth', 1)
title('Linear Regression')
ylabel('Y')
xlabel('Subjects')
legend({'Data', 'Fit'})
subplot(122)
scatter(Y, Y_hat_reg, 'r.')
xlabel('Y')
ylabel('Y_{hat}_{ridge}')
title(['R2= ' num2str(R2)])
box on

%% ESERCIZIO 2: RIDGE REGRESSION

% 6. Definire il vettore contenente i valori di coefficiente di regolarizzazione 
% da testare (usare la function logspace per definire un vettore in scala logaritmica 
% da 10-3 a 103).
lambda_ridge = logspace(-3,3,20);


% 7. Fittare i dati tramite Ridge Regression: per la scelta del parametro di 
% regolarizzazione, eseguire una Leave-One-Out-Cross-Validation, selezionando 
% il valore del coefficiente di regolarizzazione che minimizza il valore di MSE 
% (Mean Squared Error), pari alla media della somma dei residui. 

MSE = [];

for jj = 1:length(lambda_ridge) % per ogni valore jj-esimo del parametro di regolarizzazione lambda
    disp(jj)
    for ii = 1:size(data,1) % per ogni soggetto ii-esimo
        
        % Train-Test 
        X_train = X([1:ii-1 ii+1:end],:); % training: tutti i soggetti meno l'ii-esimo
        Y_train = Y([1:ii-1 ii+1:end],:);
        
        X_test = X(ii,:); % test: ii-esimo soggetto

        B_temp = (X_train'*X_train + lambda_ridge(jj).*eye(length(X_train'*X_train)))\X_train'*Y_train;
        % alternativa: function ridge
        % B_temp = ridge(Y_train,X_train,lambda(jj));
        
        Y_test_hat(ii,jj)=X_test*B_temp; % stima del soggetto ii-esimo

    end
    MSE(jj) = mean((Y-Y_test_hat(:,jj)).^2); % calcolo del MSE per il lambda jj-esimo
end

% 8. Selezionare il lambda ottimo corrispondente al minimo valore di MSE.
% Plottare l'andamento del MSE al variare del parametro di regolarizzazione 
% evidenziando il lambda ottimo. 
[~,idx_opt_ridge] = min(MSE);
lambda_OPT = lambda_ridge(idx_opt_ridge);

figure(2)
plot(lambda_ridge,MSE)
hold on
ylabel('MSE')
xlabel('parametro di regolarizzazione')
title('RIDGE REGRESSION - MSE')
plot(lambda_OPT,min(MSE),'r*')

% 9. Considerando tutti i soggetti, stimare i valori finali di beta (refit del modello) 
% usando il parametro di regolarizzazione selezionato. Calcolare la
% predizione del modello e il coefficiente di determinazione tra Y e
% predizione.
B_ridge = (X'*X+lambda_OPT.*eye(length(X'*X)))\X'*Y;
Y_hat_ridge=X*B_ridge;
R2 = corr(Y, Y_hat_ridge).^2;


% 10. Rappresentare in un plot: 
% 1) l'andamento di Y attraverso i soggetti, con sovrapposta la stima del modello 
% 2) scatter plot del dato vs la predizione.
figure(3)
subplot(121)
plot(1:n, Y, 'o--'), hold on, plot(1:n, Y_hat_ridge, 'r', 'linewidth', 1)
title(['RIDGE LOOCV, MSE =' num2str(min(MSE))])
ylabel('Y')
xlabel('Subjects')
legend({'Data', 'Fit'})
subplot(122)
scatter(Y, Y_hat_ridge, 'r.')
xlabel('Y')
ylabel('Y_{hat}_{ridge}')
title(['R2= ' num2str(R2)])
box on

% 11. [BONUS] Ripetere i punti precedenti usando 10 fold CV. Per il plot di MSE, mostrare i
% risultati di MSE sia per LOOCV, sia per 10 fold CV.
MSE_10 = [];
fold=10;
size_test=(size(data,1)/fold); % size del test set
l=size_test-1; % variabile per contare

for jj = 1:length(lambda_ridge) % per ogni valore jj-esimo del parametro di regolarizzazione lambda
    disp(jj)
    for ii = 1:50:(size(data,1)-l) % per ogni soggetto ii-esimo
        
        % Train-Test 
        X_test = X([ii:(ii+l)],:); % test size

        X_train=X;
        Y_train=Y;
        X_train([ii:(ii+l)],:) = []; % training: tutti i soggetti meno il test
        Y_train([ii:(ii+l)],:) = [];
        
        B_temp = (X_train'*X_train + lambda_ridge(jj).*eye(length(X_train'*X_train)))\X_train'*Y_train;
        % alternativa: function ridge
        % B_temp = ridge(Y_train,X_train,lambda(jj));
        
        Y_test_hat_10([ii:(ii+l)],jj)=X_test*B_temp; % stima del soggetto ii-esimo

    end
    MSE_10(jj) = mean((Y-Y_test_hat_10(:,jj)).^2); % calcolo del MSE per il lambda jj-esimo
end

% mostrare i risultati di MSE per LOOCV e 10 fold CV
[~,idx_opt_ridge_10] = min(MSE_10);
lambda_OPT_10 = lambda_ridge(idx_opt_ridge_10);

figure(4)
plot(lambda_ridge,MSE, 'b')
hold on
plot(lambda_OPT,min(MSE),'r*')
plot(lambda_ridge,MSE_10, 'g')
plot(lambda_OPT_10,min(MSE_10),'r*')
ylabel('MSE')
xlabel('parametro di regolarizzazione')
legend('LOOCV', '\lambda_{opt} LOOCV', '10=foldCV', '\lambda_{opt} 10foldCV')
title('RIDGE REGRESSION - MSE LOOCV vs 10 fold CV')

% stima dei beta, delle Y e calcolo di R2
B_ridge_10 = (X'*X+lambda_OPT_10.*eye(length(X'*X)))\X'*Y;
Y_hat_ridge_10=X*B_ridge_10;
R2_10 = corr(Y, Y_hat_ridge_10).^2;

% Rappresentare in un plot: 
% 1) l'andamento di Y attraverso i soggetti, con sovrapposta la stima del modello 
% 2) scatter plot del dato vs la predizione.

figure(5)
subplot(121)
plot(1:n, Y, 'o--'), hold on, plot(1:n, Y_hat_ridge_10, 'r', 'linewidth', 1)
title(['RIDGE 10 fold CV, MSE =' num2str(min(MSE_10))])
ylabel('Y')
xlabel('Subjects')
legend({'Data', 'Fit'})
subplot(122)
scatter(Y, Y_hat_ridge_10, 'r.')
xlabel('Y')
ylabel('Y_{hat}_{ridge}')
title(['R2= ' num2str(R2_10)])
box on

%% ESERCIZIO 3: LASSO REGRESSION

% 12. Fittare i dati tramite LASSO Regression (function lasso). 
% Settare una nuova griglia di valori per lambda (lambda=logspace(-1,0.5,200))
% Nella function usare i parametri: 'CV' = 10, 'Alpha' = 1, 'lambda' = lambda.
% Estrarre MSE e identificare il lambda ottimo. 
% Visualizzare MSE tramite lassoPlot.
lambda_vec=logspace(-1,0.5,200);

[b_lasso,fitinfo] = lasso(X,Y,'CV',10,'Alpha',1,'lambda',lambda_vec);
[~,idx_opt_lasso] = min(fitinfo.MSE);
lambda_OPT_lasso=lambda_vec(idx_opt_lasso);


% lassoPlot
figure(6)
lassoPlot(b_lasso,fitinfo,'PlotType','CV','Parent',gca);

%alternativamente, senza usare la function lassoPlot()
% figure(6)
% semilogx(lambda,fitinfo.MSE)
% hold on
% semilogx(lambda_OPT_lasso, min(fitinfo.MSE), 'r*')
% ylabel('MSE')
% xlabel('parametro di regolarizzazione (logscale)')
% legend('MSE', '\lambda_{opt}')
% title('LASSO REGRESSION - MSE')


% 13. Considerando tutti i soggetti, stimare i valori finali di beta (refit del modello) 
% usando il parametro di regolarizzazione selezionato. Calcolare la
% predizione del modello e il coefficiente di determinazione tra Y e
% predizione.
[B_lasso,fitinfoOPT] = lasso(X,Y,'CV',10,'Alpha',1,'lambda',lambda_OPT_lasso);
Y_hat_lasso=X*B_lasso;
R2_lasso = corr(Y, Y_hat_lasso).^2;

% 14. Rappresentare in un plot: 
% 1) l'andamento di Y attraverso i soggetti, con sovrapposta la stima del modello 
% 2) scatter plot del dato vs la predizione.
figure(7)
subplot(121)
plot(1:n, Y, 'o--'), hold on, plot(1:n, Y_hat_lasso, 'r', 'linewidth', 1)
title(['LASSO, MSE =' num2str(min(fitinfoOPT.MSE))])
ylabel('Y')
xlabel('Subjects')
legend({'Data', 'Fit'})
subplot(122)
scatter(Y, Y_hat_lasso, 'r.')
xlabel('Y')
ylabel('Y_{hat}_{lasso}')
title(['R2= ' num2str(R2_lasso)])
box on

%% ESERCIZIO 4: ELASTIC NET REGRESSION

% 15. Fittare i dati tramite ELASTIC NET Regression (function lasso),'CV' = 10, 'Alpha' = 0.5, 'lambda' = lambda_vec.
% Estrarre MSE e identificare il lambda ottimo. 
% Visualizzare MSE tramite lassoPlot.
[b_elnet,fitinfoElNet] = lasso(X,Y,'CV',10,'Alpha',.5,'lambda',lambda_vec);
[~,idx_opt_elastic] = min(fitinfoElNet.MSE);
lambda_OPT_elastic=lambda_vec(idx_opt_elastic);


% lassoPlot
figure(8)
lassoPlot(b_elnet,fitinfoElNet,'PlotType','CV','Parent',gca);

%alternativamente, senza usare la function lassoPlot()
% figure(8)
% semilogx(lambda,fitinfoElNet.MSE)
% hold on
% semilogx(lambda_OPT_elastic, min(fitinfoElNet.MSE), 'r*')
% ylabel('MSE')
% xlabel('parametro di regolarizzazione (logscale)')
% legend('MSE', '\lambda_{opt}')
% title('ELASTIC NET REGRESSION - MSE')

% 16. Considerando tutti i soggetti, stimare i valori finali di beta (refit del modello) 
% usando il parametro di regolarizzazione selezionato. Calcolare la
% predizione del modello e il coefficiente di determinazione tra Y e
% predizione.
[B_elastic,fitinfoElNetOPT] = lasso(X,Y,'CV',10,'Alpha',.5,'lambda',lambda_OPT_elastic);
Y_hat_elastic=X*B_elastic;
R2_elastic = corr(Y, Y_hat_elastic).^2;

% 17. Stampare nella command window i valori minimi di MSE ottemito da lasso e Elastic Net 
disp(['Il valore minimo di MSE ottenuto con lasso è: ' num2str(min(fitinfo.MSE))])
disp(['Il valore minimo di MSE ottenuto con Elastic NET è: ' num2str(min(fitinfoElNet.MSE))])

% codice per il plot manuale dell'andamento di MSE per lasso ed Elastic NET
% figure(9)
% semilogx(lambda_vec,fitinfo.MSE)
% hold on
% semilogx(lambda_OPT_lasso, min(fitinfo.MSE), 'r*')
% semilogx(lambda_vec,fitinfoElNet.MSE)
% semilogx(lambda_OPT_elastic, min(fitinfoElNet.MSE), 'r*')
% ylabel('MSE')
% xlabel('parametro di regolarizzazione (logscale)')
% legend('MSE LASSO', '\lambda_{opt} LASSO', 'MSE ELASTIC', '\lambda_{opt} ELASTIC')
% title('LASSO VS ELASTIC NET - MSE')

% 18. Rappresentare in un plot: 
% 1) l'andamento di Y attraverso i soggetti, con sovrapposta la stima del modello 
% 2) scatter plot del dato vs la predizione.
figure(9)
subplot(121)
plot(1:n, Y, 'o--'), hold on, plot(1:n, Y_hat_elastic, 'r', 'linewidth', 1)
title(['ELASTIC, MSE =' num2str(min(fitinfoElNet.MSE))])
ylabel('Y')
xlabel('Subjects')
legend({'Data', 'Fit'})
subplot(122)
scatter(Y, Y_hat_elastic, 'r.')
xlabel('Y')
ylabel('Y_{hat}_{elastic}')
title(['R2= ' num2str(R2_elastic)])
box on

%% ESERCIZIO 5: CONFRONTO
% 19. Plottare i Beta stimati tramite ridge regression LOOCV, quelli stimati
% tramite lasso regression e quelli dell'elastic net. Cosa si può dire rispetto alla selezione delle
% variabili da includere nel modello? 
figure(10)
stem(1:1:size(X,2), abs(B_reg),'k*')
hold on
stem(1:1:size(X,2), abs(B_ridge), 'b*')
hold on
stem(1:1:size(X,2), abs(B_lasso), 'r*')
hold on
stem(1:1:size(X,2), abs(B_elastic), 'g*')
legend('\beta REG','\beta RIDGE', '\beta LASSO', '\beta ELASTIC NET')
xlabel('variables')
ylabel('\beta')
title('\beta Plot')


disp('# of zero-coefficients with ridge-regression:')
disp(sum(B_ridge==0))
disp('# of zero-coefficients with lasso:')
disp(sum(B_lasso==0))
disp('# of zero-coefficients with elastic-net:')
disp(sum(B_elastic==0))