%%
clear all
close all
clc
%%
load dati_Ridge


%%
% stima ai minimi quadrati con varianza non nota da stimare

stime=inv(X'*X)*X'*Y;
predizione=X*stime;

% valutazione della bontà del fit tramire regressione
[rho,pval] = corr(Y,predizione);
R2=rho*rho;
res=Y-predizione;

sigma2=(res'*res)/(length(Y)-length(stime));
cov_err=(sigma2).*inv(X'*X);
figure, subplot(221), plot(res), title('RESIDUI')
sd_perc_err=sqrt(diag(cov_err))./abs(stime)*100;
subplot(223), bar(sd_perc_err), ylim([0 100]), title('errore di stima %')
subplot(222), plot(predizione,'k-','LineWidth',1), hold on, plot(Y,'og','LineWidth',2), xlim([0 204])
subplot(224), bar(stime)

clc
R2=corr(Y,X*stime)^2;
disp(['R2 = ', num2str(R2)])

%%
%lambda=0.1, 10, 100000
lambda=10000;
B = ridge(Y,X,lambda);
stime_R=B;

figure, stem(B,'r'), hold on, stem(stime,'k')

predizione_R=X*B;

figure
subplot(311), plot(predizione_R,'-r'), hold on, plot(Y,'bo'),title('confronto predizioni Ridge'), xlim([0 205]),xlabel('Ridge'), ylabel('minimi quadrati')
subplot(312), plot(predizione,'-k'), hold on, plot(Y,'bo'),title('confronto predizioni LS'), xlim([0 205]),xlabel('Ridge'), ylabel('minimi quadrati')
subplot(313), stem(stime_R,'r'), hold on, stem(stime,'k'),  title('stime'), ylim([-1 1]) %ylim([min(stime) max(stime)])
%%
best_lambda = ridge_regression_lambda_loocv(X, Y)

B_best_lambda = ridge(Y,X,best_lambda);
stime_R_best_lambda=B_best_lambda;

figure, stem(B_best_lambda,'r'), hold on, stem(stime,'k')

predizione_R_best_lambda=X*B_best_lambda;

figure
subplot(311), plot(predizione_R_best_lambda,'-r'), hold on, plot(Y,'bo'),title('confronto predizioni Ridge'), xlim([0 205]),xlabel('Ridge'), ylabel('minimi quadrati')
subplot(312), plot(predizione,'-k'), hold on, plot(Y,'bo'),title('confronto predizioni LS'), xlim([0 205]),xlabel('Ridge'), ylabel('minimi quadrati')
subplot(313), stem(stime_R_best_lambda,'r'), hold on, stem(stime,'k'),  title('stime'), ylim([min(stime) max(stime)])
%%
% Esempio di stima elastic-net 
lambda=1;
B = lasso(X,Y,'Alpha',0.2,'Lambda',lambda,'Intercept',false)
predizione=X*B;
R2=corr(Y,predizione)^2

figure, plot(Y,'k'), hold on, plot(predizione,'-m')