clear
close all
clc

% Load del dato
load ('data_SUVR.mat')

%% Pulizia del dato
flag_2_delete_G = any(isnan(G_bivar) | isinf(G_bivar) | (G_bivar<0),2);
G_bivar(flag_2_delete_G,:) = [];
SUVR(flag_2_delete_G) = [];

% N.B. controllo anche la variabile dipendente SUVR ed elimino i soggetti
% corrispondenti
flag_2_delete_SUVR = any(isnan(SUVR) | isinf(SUVR) | (SUVR<0),2);
G_bivar(flag_2_delete_SUVR,:) = [];
SUVR(flag_2_delete_SUVR) = [];

%% Z-score
Y=zscore(SUVR);
X=zscore(G_bivar);

%% Pre-alloco le variabili per essere piu' efficente 
rss      = zeros(size(X,2)-1,1);
res      = zeros(size(Y,1),size(X,2)-1);
y_hat    = zeros(size(Y,1),size(X,2)-1);
rho      = zeros(size(X,2)-1,1);
R2       = zeros(size(X,2)-1,1);
p_val_R2 = zeros(size(X,2)-1,1);

% G = sottoinsieme delle colonne di X:
 % i=2 -> elimino le colonne da 2 in poi
 % i=end -> elimino solo l'ultima colonna
for i=2:size(X,2)+1 %il +1 ci permette di considerare poi anche l'ultima colonna di X
    
    % Per ogni G preparo i calcoli dei parametri e i grafici da mostrare
    G=X;
    % Elimino le colonne da i in poi. Se i=size(X,2) prendo tutta a matrice
    % X quindi non devo eliminare nessuna colonna
    if i<size(X,2)
        G(:,i:end)=[];
    end
   
    % Calcolo i parametri beta con il modello aggiornato che comprende le
    % i-1 variabili (colonne di X)
    beta_hat = (G'*G)\G'*Y;
    
    % Calcolo la predizione del modello 
    y_hat(:,i-1) = G*beta_hat;

    % Calcolo i residui
    res(:,i-1) = Y-y_hat(:,i-1);
    
    % Calcolo la somma dei residui al quadrato
    rss(i-1) = (res(:,i-1)'*res(:,i-1));

    figure(1)
    clf
    subplot(221)
    plot(Y,'ko')
    hold on
    plot(y_hat(:,i-1))
    xlim([1 size(res,1)])
    legend('Dati','Predizione')
    title('Fit')

    subplot(222)
    plot(res(:,i-1))
    xlim([1 size(res,1)])
    ylim([-2 2])
    title('Residui')

    subplot(223)
    plot(rss)
    hold on
    plot(i-1,rss(i-1),'r+')
    xlim([1 52])
    title('Root Sum of Squares')

   
    sigma_hat_2 = (res(:,i-1)'*res(:,i-1))/(length(Y)-length(beta_hat));
    
    % Calcolo lo standard error
    SE =sqrt((sigma_hat_2)*diag(inv(G'*G))); 

    % Calcolo i coefficienti di variazione
    CV = SE./abs(beta_hat)*100;

    subplot(224)
    stem(CV)
    xlabel('Variabili indipendenti')
    ylabel('CV[%]')
    xlim([1 52])
    ylim([0 100])
    title('Errore delle stime dei beta')

    % Calcolo il coeff. di determinazione R2
    % i.e. per modelli lineari e' equivalente alla correlazione al quadrato
    [rho(i-1),p_val_R2(i-1)] = corr(Y,y_hat(:,i-1));
    % N.B. il secondo parametro e' il p-value, che indica se la
    % correlazione e' significativa
    R2(i-1)=rho(i-1).^2;

    drawnow()
    pause(0.5)
    % calcolo la stima a posteriori dell'errore di misura (N.B. cambia in
    % relazione ai residui e al numero di beta, i.e. numero di parametri
    % del modello)
    
end

%% Risultati modello con una variabile
figure(2)
subplot(231)
plot(Y,'ko')
hold on
plot(y_hat(:,1))
xlim([1 size(res,1)])
legend('Dati','Predizione')
title('Fit')

subplot(232)
scatter(Y,y_hat(:,1))
xlabel('Dati (Y)')
ylabel('Predizione (Yhat)')
title('Dati vs Predizione')

subplot(233)
plot(res(:,1))
title('Residui')
xlim([1 size(res,1)])

subplot(234)
autocorr(res(:,1),'NumSTD',2)  

subplot(235)
qqplot(res(:,1))

%% Risultati modello con matrice X completa
figure(3)
subplot(231)
plot(Y,'ko')
hold on
plot(y_hat(:,end))
xlim([1 size(res,1)])
legend('Dati','Predizione')
title('Fit')

subplot(232)
scatter(Y,y_hat(:,end))
xlabel('Dati (Y)')
ylabel('Predizione (Yhat)')
title('Dati vs Predizione')

subplot(233)
plot(res(:,end))
title('Residui')
xlim([1 size(res,1)])

subplot(234)
autocorr(res(:,end),'NumSTD',2)  

subplot(235)
qqplot(res(:,end))


%% Plot di Rss al variare di G
%  parametri beta e CV del modello
figure(4)
subplot(311)
plot(rss)
xlim([1 52])
title('Somma dei residui al quadrato')
subplot(312)
stem(beta_hat)
xlim([1 52])
title('Coeff. Beta')
subplot(313)
stem(CV)
title('Coeff. di variazione percentuali')
xlim([1 52])
ylim([0 100])

%% Calcolo di R2 e relativo p-value per il modello a una variabile e quello con la matrice X completa

disp(['R2 - X con singola variabile:' num2str(R2(1))])
disp(['p-value ' num2str(p_val_R2(1))])

disp(['R2 - X con tutte le variabili:' num2str(R2(end))])
disp(['p-value ' num2str(p_val_R2(end))])