%% Testo
% 1) Caricare  i dati presenti nel file di Excel conadR.xlsx
% 2) Calcolare e commentare la matrice di correlazione ed i relativi pvalue del test H_0 :\rho =0\text{ }
% 3) Calcolare manualmente la matrice dei p-values
% 4) Verificare empiricamente la distribuzione T di student con n-2 gradi di libertà del test sul coefficiente di correlazione lineare confrontando i quantili empirici e quelli teorici.
% Utilizzare i quantili che seguono.
% quant=[0.01 0.05 0.10:0.1:0.9 0.95 0.99 0.999];
% 5) Rappresentare graficamente la soglia oltre il quale il valore di rxy è ritenuto significativo al livello del 5% in funzione della numerosità campionaria
% 6) Caricare i dati presenti nel file tagliatelle.xlsx. Importare i dati in una matrice denominata X.
% 7) Rimuove in maniera listwise i valori mancanti. Sulla matrice risultatante calcolare le correlazioni ed i reltivi p-values
% 8) Rimuovere dalla matrice X  le colonne che hanno almeno 5 missing values.
% 9) Utilizzando la funzione nancov sulla matrice X calcolare la matrice di covarianza. Partendo dalla matrice di covarianza calcolare manualmente la matrice di correlazione.
% 10) Utilizzando l'opzione 'rows' della funzione corr confrontare la matrice di correlazione ottenuta tramite cancellazione listwise e cancellazione pairwise dei dati mancanti.
% 11) Caricare dentro MATLAB i dati relativi ai passeggeri e merci di diversi aeroporti. Calcolare e commentare il coefficiente di cograduazione di Spearman. Verificare il risultato ottenuto passando attraverso la funzione tiedranks. Verificare il risultato implementando direttamente la formula 1- 6 \sum .....
%%
%% Soluzione.
%% Carico i dati del file conadR.xlsx dentro MATLAB

[X] = xlsread('conadR.xlsx','dati','A1:D509');

% Calcolo della matrice di correlazione e dei p-values del test
[R,Pval]=corr(X);
disp('Matrice di correlazione')
disp(R)
disp('Matrice dei p-values')
disp(Pval)

%
%% Calcolo manuale della matrice dei p-values

% Calcolo del test per ogni coppia di variabili
Testt=(R./sqrt(1-R.^2))*sqrt(n-2);

% Calcolo dei relativi p-values
% Moltiplico l'area della coda di destra per 2 poiché l'ipotesi alternativa
% è bilaterale
Pvalchk=tcdf(abs(Testt),n-2,'upper')*2;

%
%% Verifica empirica della distribuzione T di student con n-2 gdl

n=10;
nsimul=100000;
Test=zeros(nsimul,1);
for i=1:nsimul
    x=randn(n,2);
    [R]=corr(x);
    R=R(1,2);
    % Inserisci il valore del test nella riga i del vettore Test
    Test(i)=(R/sqrt(1-R^2))*sqrt(n-2);
end
% Valori ordinati
Testsor=sort(Test);
quant=[0.01 0.05 0.10:0.1:0.9 0.95 0.99 0.999]';
Empirici=Testsor(round(nsimul*quant));
Teorici=tinv(quant,n-2);
disp('Confronto tra quantili empirici e quantili teorici')
disp([quant Empirici Teorici]);

% Confronto tra quantili empirici e quantili teorici
plot(Empirici,Teorici,'o')
xlabel('Quantili empirici')
ylabel('Quantili teorici')

% quantili ottenuti tramite la funzione prctile
% plot(prctile(Test,100*quant),Teorici,'o')

%
%% Valore soglia per la significatività di r al variare di n
% Rappresentare graficamente la soglia oltre il quale il valore di rxy è ritenuto significativo al livello del 5% in funzione della numerosità campionaria

nn=5:100;
tinvalpha2=(tinv(0.975,nn-2)).^2;
soglia=sqrt(tinvalpha2./(nn-2+tinvalpha2));
plot(nn,soglia)
ylabel('Valore soglia per la significatività di r al 5%')
xlabel('Numerosità campionaria')

%
%% TRATTAMENTO DEI DATI MANCANTI

% Carico i dati dentro MATLAB
[X,nomi] = xlsread('tagliatelle.xlsx','dati','A1:C41');
% rimozione listwise dei valori mancanti
Xnomiss=rmmissing(X,1);
[n,p]=size(Xnomiss);
disp('Numero di righe senza missing values')
disp(n)
% Calcolo della matrice di correlazione e dei relativi p-values
[R,Pval]=corr(Xnomiss);
disp('Matrice di correlazione')
disp(R)
disp('Matrice dei p-values')
disp(Pval)

%
%% Rimuovere le colonne che hanno un certo numero di missing
% Rimuovi le colonne che hanno almeno 5 missing values

[Rred,TF] = rmmissing(X,2,'MinNumMissing',5)

%
%% Matrice dei coefficienti di correlazione usando nancov e un ciclo
% Calcolo la matrice di covarianze escludendo le righe con almeno un valore uguale a NaN

covY = nancov(X);

% R = matrice pxp di correlazione
Rchk=zeros(p,p);
for i=1:p
    for j=1:p
        Rchk(i,j)=covY(i,j)/sqrt(covY(i,i)*covY(j,j));
    end
end

%
%% Calcolo della matrice di correlazione utilizzando l'opzione dentro corr

[Rlistwise,Pvallistwise]=corr(X,'rows','complete');

[Rpaiwise,Pvalpairwise]=corr(X,'rows','pairwise');
disp('Confronto tra le due matrici di correlazione')
disp('Matrice di correlazione dopo dopo l''eliminazione listwise')
disp(Rlistwise)
disp('Matrice di correlazione dopo dopo l''eliminazione pairwise')
disp(Rpaiwise)

%
%% Calcolo della cograduazione

eti={'Bergamo-Orio'
    'Bologna-BorgoP'
    'Cagliari-Elmas'
    'Milano-Linate'
    'Milano-Malpensa'
    'Roma-Fiumicino'
    'Venezia-Tessera'
    'Verona-Villafranca'};
dati=  [ 47820    4291239;
    54780   3624072;
    26425   2344282;
    93942   9085999;
    227718  19499158;
    302890  28208161;
    75196   5780783;
    33178   2581420];

[R,Pval]=corr(dati,'type','Spearman');
disp('Matrice di cograduzione')
disp(R)

%
%% Implementazione utilizzando la funzione tiedrank

datirank=zeros(size(dati));
for j=1:2
    datirank(:,j)=tiedrank(dati(:,j));
end
Rchk=corr(datirank)

%
%% Implementazione utilizzando direttamente la formula in funzione dei ranks

n=size(datirank,1);
R12chk1=1-6*sum((datirank(:,1)-datirank(:,2)).^2)/(n*(n^2-1));
% Controllo sulle due implementazioni
disp(Rchk(1,2)-R12chk1)

%
%%