% Verificare empiricamente la distribuzione T di student con n-2 gradi 
% di libertà del test sul coefficiente di correlazione lineare confrontando i 
% quantili empirici con quelli teorici. Utilizzare 100000 simulazioni ed un
% valore di n pari a 10 
% 
% Utilizzare i quantili che seguono.
% 
% quant=[0.001	0.01	0.05	0.1	0.15	0.2	0.25	0.3	0.35	0.4 	0.45	0.5	0.55	0.6	0.65	0.7	0.75	0.8	0.85	0.9	0.95	0.99	0.999];
% 
n=10;
nsimul=100000;
% Test = vettore di dimensione 100000x1 che conterrà in posizione i il test
% sul coefficiente di correlazione basato sulla simulatione i
Test=zeros(nsimul,1);
% Rall= vettore di dimensione nsimul x1 che contiene in posizione i il
% valore del coefficiente di correlazione per la simulazione i-esima
Rall=zeros(nsimul,1);
for i=1:nsimul
    % Per ogni simulazione vengono generate due variabili indipendenti di
    % dimensione nx1
    X=randn(n,2);
    % Per ogni simulazione si calcola il coefficiente di correlazione
    % lineare
    [R]=corr(X);
    R=R(1,2);
    % Inserisco il valore del test per la simulazione i nella riga i del
    % vettore Test
    Test(i)=(R/sqrt(1-R^2))*sqrt(n-2);
    % Inserisco il valore del coefficiente di correlazione per la
    % simulazione i nella riga i del vettore Rall
    Rall(i)=R;
end
% Valori ordinati
Testsor=sort(Test);
quant=[0.001	0.01	0.05	0.1	0.15	0.2	0.25	0.3	0.35	0.4 	0.45	0.5	0.55	0.6	0.65	0.7	0.75	0.8	0.85	0.9	0.95	0.99	0.999]';

%% Calcolo dei quantili empirici con quelli teorici dell T(n-2)
Empirici=Testsor(round(nsimul*quant));
% Calcolo dei quantili teorici utilizzando la T di student con (n-2) gradi
% di libertà
Teorici=tinv(quant,n-2);
disp('Confronto tra quantili empirici e quantili teorici')
disp([quant Empirici Teorici]);

% Costruire un grafico a dispersione che contiene sull'asse delle ascisse i
% quantili empirici e sull'asse delle ordinate i quantili teorici
%
% Confronto tra quantili empirici e quantili teorici in maniera grafica
plot(Empirici,Teorici,'o')
xlabel('Quantili empirici')
ylabel('Quantili teorici della T(n-2)')


%% P value empirici e teorici della distribuzione N(0,1)
TeoriciN=norminv(quant);
% Confronto tra quantili empirici e quantili teorici in maniera grafica
plot(Empirici,TeoriciN)
xlabel('Quantili empirici')
ylabel('Quantili teorici N(0,1)')

% La distribuzione N(0,1) presenta una probabilità molto più bassa nelle
% code della T(8) di conseguenza il grafico che riporta sull'asse delle
% ordinate i percentili teorici della N(0,1) e sull'asse delle ascisse i
% percentili della normale standardizzata presenta concavità rivolta verso
% l'alto nell'intervallo [-5 0] in quanto, in questo intervallo i decili
% della T(n-2) sono più piccoli di quelli della N(0,1) e concavità rivolta
% verso il  basso nell'intervallo [0 5] in quanto in questo intervallo i
% decili della T(8) sono più grandi di quelli della N(0,1)

%% 
% Supponiamo di avere estratto un campione casuale bivariato di 10 elementi
% e che il test sulla correlazione (H_0:\rho=0) sia risultato pari a -2.6.
% Calcolare il p value empirico e teorico del numero -2.6 in presenza di:
% 1) ipotesi alternativa bilaterale H_1: \rho <> 0 (\rho diverso da zero)
pvalEmpiricoH1bil=sum(abs(Test)>2.6)/nsimul;
pvalTeoricoH1bil=tcdf(abs(-2.6),n-2,'upper')*2;
disp('p-value teorico calcolato tramite la distribuzione T(n-2) H_1: \rho <> 0 ')
disp(pvalTeoricoH1bil)
disp('p-value teorico calcolato estraendo 100000 campioni')
disp(pvalEmpiricoH1bil)

% 2) ipotesi alternativa unilaterale destra H_1: \rho > 0
pvalEmpiricoH1dx=sum(Test>-2.6)/nsimul;
pvalTeoricoH1dx=tcdf(-2.6,n-2,'upper');
disp('p-value teorico calcolato tramite la distribuzione T(n-2) H_1: \rho <> 0 ')
disp(pvalTeoricoH1dx)
disp('p-value teorico calcolato estraendo 100000 campioni')
disp(pvalEmpiricoH1dx)

% 3) ipotesi alternativa unilaterale sinistra. H_1: \rho < 0
pvalEmpiricoH1sx=sum(Test<-2.6)/nsimul;
pvalTeoricoH1sx=tcdf(-2.6,n-2);
disp('p-value teorico calcolato tramite la distribuzione T(n-2) H_1: \rho <> 0 ')
disp(pvalTeoricoH1sx)
disp('p-value teorico calcolato estraendo 100000 campioni')
disp(pvalEmpiricoH1sx)

% Osservazione finale: nella pratica io non osservo i risultati
% dell'applicazione del test su 100000 campioni, in quanto   il risultato
% del test è basato su un solo campione. Dato che il risultato del test si
% distribuisce come una v.c. T di Student con (n-2) gradi di libertà, io
% posso sapere utilizzando la funzione di ripartizione della v.c. T di
% Student con (n-2) gradi di libertà la probabilità del risultato
% campionario ossia la probabilità che l'ipotesi nulla sia vera. In
% presenza di ipotesi alternativa bilaterale la probabilità del risultato
% campionario è la probabilità di ottenere un valore in modulo superiore a
% quello del test ... Probabilità bassa del risultato campionario (ossia
% p-value basso) mi portano a ritenere che l'ipotesi nulla di assenza di
% correlazione non sia accettibile. Quindi, in presenza di p-value basso
% accetto l'ipotesi alternativa di correlazione diversa da zero tra le due
% variabili in esame.