%% PARTE I
% Fissare il seed dei numeri casuali a 100.
% Assegnare ad un variabile n il valore 200 (n è la variabile che contiene
% la numerosità campionaria)
% Generare i dati x in base alla distribuzione uniforme nell'intervallo [-5 5] e i disturbi epsilon in base alla distribuzione normale N(0, sigma).
% Fissare alpha=12, beta=3 e sigma=10
%
% Stimare i parametri a e b del modello di regressione e l'R2 tramite la funzione di MATLAB fitlm.
% Verificare il risultato ottenuto per a e b tramite l'implementazione manuale
% utilizzando per la stima dell'intercetta e della pendenza le formule
% riportate di seguito.



%
%
%
% Calcolare la stima dei parametri a e b della retta di regressione
% direttamente tramite un vettore di dimensione 2x1 utilizzando
% utilizzando bvec=X\y; Costruire in maniera appropriata la matrice X
%
% Costruire il diagramma di dispersione con sovrapposta  la retta di regressione.
% Inserire nel titolo del grafico il valore dei parametri stimati e di R2.



%
%
%% PARTE II



% Esplorare cosa succede se si cambiano i parametri alpha, beta e sigma.
% Copiare il codice utilizzato sopra per generare i dati.
% Contaminare la prima osservazione aggiungendo un valore definito nella variabile
% valorecont.
% Associare ai parametri alpha, beta, n, sigma e valorecont 5 sliders
% Lo slider di sigma deve avere come valore min 0.0001, max 20 e step 0.1
% Lo slider di n deve avere come valore min 20 max 400 e step 10
% Lo slider di alpha deve avere come valore min 10, max 30 e step 1
% Lo slider di beta deve avere come valore min 0, max 10 e step 1
% Lo slider di valorecont deve avere come min 0 max 2000 e step 1
%
% Il teorema di Sastry e Nag afferma che la presenza di un outlier (osservazione molto distante dalla nuvola dei punti) fa andare ad 1 il valore di R2. Verificare empiricamente cosa succede alla stima di R2 quando valorecont aumenta e tende all'infinito.
% Esplorare l'opzione "Hide code" nella barra laterale destra.



%
%% Soluzione generazione di x, eps e y

% fisso il seed dei numeri casuali
rng(100)
% assegno alla variabile n il valore 200
n=200;
% Fisso i parametri richiesti
sigma=10; % dispersione dei valori attorno alla retta di regressione
alpha=12; % valore della vera intercetta
beta=3; % valore della vera pendenza della retta di regressione

% vettore x (variabile explicativa simulata) generato in base ad una
% distribuzione U(-5, 5)
x=unifrnd(-5,5,n,1);
% disturbi espilon generati in base a N(0, sigma)
eps=normrnd(0,sigma,n,1);
% Genero le osservazioni y in base al modello y=alpha+beta*x+eps
y=alpha+beta*x+eps;

%
%% Stima dei parametri del modello di regressione tramite la funzione fitlm

out=fitlm(x,y);
disp(out)

% Osservazione: class(out) produce 'LinearModel'
% out è un classe particolare di tipo "linearModel che contiene determinate
% proprietà e metodi (funzioni)
% Ad esempio out.Coefficients è una table con 4 colonne che contiene le
% stime dei parametri, gli standard error, le statistiche t ed i relativi
% p-value
disp(out.Coefficients)
% Estraggo dalla table out.Coefficients i valori di alpha cappello e beta
% cappello ossia dell'intercetta e della pendenza stimati con il metodo dei
% minimi quadrati ed il valore di R2 (coefficiente di determinazione).
alphastimato=out.Coefficients{1,1};
betastimato=out.Coefficients{2,1};
R2=out.Rsquared.Ordinary;
% Grafico dei punti x y 
plot(x,y,'o')
xlabel('Dati x simulati in base a U(-5,5)')
% Il comando lsline sovrappone ai punti la retta di regressione
h=lsline;
h.Color='r';
title(['Parametri stimati a=' num2str(alphastimato) ...
    ' b=' num2str(betastimato) ' R2=' num2str(R2)])

%% Implementazione di a e b utilizzando le formule
%
%

sumy=sum(y);
sumx=sum(x);
sumxy=sum(x.*y);
sumx2=sum(x.^2);
sumxsq=(sum(x))^2;

denom=n*sumx2-sumxsq;
% numa = numeratore di a
numa=sumy*sumx2-sumx*sumxy;
% numb=numeratore di b
numb=n*sumxy-sumx*sumy;

%

a=numa/denom;
b=numb/denom;
disp('Stima intercetta implementazione manuale')
disp(['a=' num2str(a)])
disp('Stima pendenza implementazione manuale')
disp(['b=' num2str(b)])

%
%%  Calcolo di a e b direttamente tramite l'istruzione bvec=X\y;

% In maniera matriciale il modello y_i=a+bx_i + epsilon_i
% si scrive y=X*bvec+epsilon
% dove y= nx1
% X= nx2
% epsilon= nx1


X=[ones(n,1) x];
% Il primo elemento di bvec contiene a (stima dell'intercetta) mentre il secondo elemento di di
% bvec contiene b (stima della pendenza) 
bvec=X\y;
disp('Stime dei coefficienti calcolati in maniera matriciale')
disp(bvec)

% L'obiettivo della seconda parte è quello di capire cosa succede alla
% retta di regressione quando cambio i parametri iniziali (n, sigma, alpha, beta) con cui ho
% generato le diverse quantità


%
%% Come cambia la stima della retta di regressione quando cambio i parametri iniziali


n=120;
% Fisso i parametri
sigma=3; 
alpha=10; % valore della vera intercetta
beta=1; % valore della vera pendenza della retta di regressione

valorecont=0; % valore di contaminazione per la prima ossservazione

x=unifrnd(-5,5,n,1); eps=normrnd(0,sigma,n,1);
y=alpha+beta*x+eps; 
x(1)=x(1)+valorecont;
y(1)=y(1)+valorecont;
out=fitlm(x,y);
alphastimato=out.Coefficients{1,1};
betastimato=out.Coefficients{2,1};
R2=out.Rsquared.Ordinary;
% Grafico dei punti x y 
plot(x,y,'o')
xlabel('Dati x simulati in base a U(-5,5)')
h=lsline;
h.Color='r';
title(['Parametri stimati a=' num2str(alphastimato) ...
    ' b=' num2str(betastimato) ' R2=' num2str(R2)])

% Commento: all'aumentare del valore di sigma la variabilità degli errori aumenta e l'R2 stimato diventa sempre più piccolo
% La presenza di un outlier (osservazione molto distante dalla nuvola dei punti) fa andare ad 1 il valore di R2 (verifica empirica del teorema di Sastry e Nag).
%



