%% Caricare in memoria i dati sull'andamento della popolazione in USA
% (dataset census).
% Mostrare il grafico dell'andamento della popolazione
%
% Adattare una funzione interpolante lineare tramite la funzione fitlm.
% Costruire il grafico dei valori effettivi, dei valori adattati e di un
% intervallo di confidenza dei valori adattati.
%
% Analizzare il grafico dei residui di Pearson
% definiti come (y-yhat)/sigmahat tramite la funzione resindexplot di FSDA
% toolbox dove sigmahat è la stima non distorta dello scostamento
% quadratico medio degli errori. Commentare l'andamento dei residui al
% variare del tempo
%
% Verificare che la stima dei parametri è
% esattamente quella che si ottiene
% tramite l'operatore \ costruendo
% manualmente la matrice delle variabili esplicative.
%
% Ripetere tutti i passaggi precedenti adattando un trend di tipo
% quadratico.
% Scrivere un loop in modo da rappresentare graficamente l'andamento della
% bontà di adattamento all'aumentare del grado del polinomio interpolante.
% Fermarsi al terzo grado. Commentare il grafico che si ottiene.
%


%% Carico i dati richiesti
load census
y=pop;
plot(cdate,y)
xlabel('Anni')
ylabel('Andamento della popolazione')

%% Adattamento di una funzione interpolante lineare tramite la funzione fitlm
X=zscore(cdate);
out=fitlm(X,pop);
disp(out)

%% E' possibile chiamare la funzione plot direttamente sull'oggetto prodotto da fitlm
% L'istruzione plot chiamata sull'output di fitlm consente di ottenere
% direttamente il grafico dei valori effettivi, dei valori adattati e di un
% intervallo di confidenza dei valori adattati
plot(out)

%% Stima dei parametri della funzione lineare tramite operatore \
% Costruzione della matrice X con intercetta
Xwithintercept=[ones(length(cdate),1) X];
b=Xwithintercept\pop;
disp('Parametri stimati tramite l''operatore \')
disp(b)

%% Analisi dei residui di Pearson
res=out.Residuals{:,2};
% oppure
% res=out.Residuals.Pearson;
resindexplot(res)


%% Adattamento di una funzione quadratica tramite la funzione fitlm
% MODO 1: si costruisce direttamente la matrice delle variabili esplicative
x1=zscore(cdate);
X=[x1 x1.^2];
out=fitlm(X ,pop);
disp(out)

%% Adattamento di una funzione quadratica tramite la funzione fitlm
% MODO 2: si chiama fitlm con l'opzione quadratic
out1=fitlm(x1 ,pop,"quadratic");
disp(out1)

%% Intervallo di confidenza dei valori adattati
plot(out1)

%% Plot dei residui utilizzando un trend quadratico
resindexplot(out1.Residuals{:,2},'conflev',0.99)

%% Stima dei parametri della funzione quadratica tramite operatore \
% Costruzione della matrice X con intercetta
Xwithintercept=[ones(length(cdate),1) X];
b=Xwithintercept\pop;
disp('Parametri stimati tramite l''operatore \')
disp(b)


%% Scelta del grado del miglior polinomio interpolante
close all
seqpol=1:3;

Gof=[seqpol' zeros(length(seqpol),1)];
X=[x1 x1.^2 x1.^3];


for i=seqpol
   out=fitlm(X(:,1:i),pop);
    Gof(i,2)=out.Rsquared.Ordinary;
end

plot(Gof(:,1),Gof(:,2),'--r')
xlabel('Ordine del polinomio interpolante')
ylabel('$R^2$','Interpreter','latex');
title('Analisi del valore di $R^2$ al variare del grado di adattamento',...
    'Interpreter','latex')
set(gca,'XTick',Gof(:,1))