%% Testo
% 1) Caricare  i dati presenti nel file di Excel conadR.xlsx sia con la
% funzione xlsread sia tramite readtable. Controllare che i dati numerici
% contenuti dentro il file di Excel vengano caricati nello stesso modo
% nelle due implementazioni.
% Costruire la matrice dei diagrammi di dispersione delle 4 variabili in
% esame prima utilizzando la funzione gplotmatrix e poi la funzione spmplot
% di FSDA toolbox. Inserire i boxplot lungo la diagonale principale nella
% matrice dei diagrammi di dispersione. Commentare la scatter plot matrix.
% 
% 2) Calcolare e commentare la matrice di correlazione ed i relativi pvalue
% del test $H_0 :\rho =0$. Mostrare la matrice di correlazione e la matrice
% dei p-values in formato table.
% 
% 3) Calcolare manualmente la matrice dei p-values.
%    Interpretare il p-value tra le variabili "visite" ed "età".
% 
% 4) In un campione bivariato di 12 pezzi da un processo produttivo il
%   coefficiente di correlazione è risultato pari a 0.54.
%   Verificare che il p-value teorico del test (nell'ipotesi di ipotesi
%   alternativa bilateale) è pari a 0.0699.
%   Se invece di un solo campione di 12 pezzi noi disponessimo di 100000
%   campioni di 12 pezzi, se è vera l'ipotesi nulla di assenza di
%   correlazione tra le due variabili su quanti campioni ci aspettiamo di un
%   ottenere un valore del coefficiente di correlazione superiore a 0.54 in
%   modulo?
%
% 5) 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 100000
%   simulazioni ed un valore di n a piacere.
% 
% Utilizzare i quantili che seguono.
% 
% quant=[0.01 0.05 0.10:0.1:0.9 0.95 0.99 0.999];
% 
% Costruire un grafico a dispersione che contiene sull'asse delle ascisse i
% quantili empirici e sull'asse delle ordinate i quantili teorici.
%
% 6) Nel caso di n=12  sulle 100000 repliche
% quanti sono le volte in cui il coefficiente di correlazione risulta
% superiore a 0.54? Confrontare la frequenza relativa empirica del numero
% di volte in cui |r|>0.54 con il p-value teorico calcolato al punto 4)
%
% 7) Rappresentare graficamente la soglia oltre il quale il valore di rxy è
% ritenuto significativo al livello del 5% in funzione della numerosità
% campionaria

%% Soluzione.

%% Carico i dati del file conadR.xlsx dentro MATLAB
% 1) Caricare  i dati presenti nel file di Excel conadR.xlsx.
[X,nomi] = xlsread('conadR.xlsx','dati','A1:D509');
% Modo alternativo per caricare i dati
Xtable=readtable('CONADR.xlsx','VariableNamingRule','preserve');

% Controllo che i due modi di caricare i dati producano gli stessi
% risultati

if isequal(Xtable{:,:},X)
else
    disp('Attenzione le istruzioni xlsread e readtable producono risultati diversi')
end

%% Calcolo della matrice dei diagrammi di dispersione con gplotmatrix
gplotmatrix(X,[],[],[],[],[],[],[],nomi);

%% Calcolo della matrice dei diagrammi di dispersione tramite spmplot
spmplot(Xtable);

%% Calcolo della matrice dei diagrammi di dispersione tramite spmplot
spmplot(Xtable,'dispopt','box');

%% Calcolo della matrice di correlazione e dei p-values del test
% 2) Calcolare e commentare la matrice di correlazione ed i relativi pvalue
% del test $H_0 :\rho =0$. Mostrare la matrice di correlazione e la matrice
% dei p-values in formato table.

[R,Pval]=corr(X);
disp('Matrice   di correlazione')
disp(R)
array2table(R,'RowNames',nomi,'VariableNames',nomi)


disp('Matrice dei p-values')
disp(Pval)
% Matrice dei coefficienti di correlazione in formato table
disp('Matrice dei coefficienti di correlazione')
disp(array2table(R,'VariableNames',nomi,'RowNames',nomi))

disp('Matrice dei pvalues')
disp(array2table(Pval,'VariableNames',nomi,'RowNames',nomi))

% Commento: al livello di significatività del 5% tutte le correlazioni sono
% significative (i p values sono tutti inferiori a 0.05).
% Da segnalare che il pvalue relativo all'assenza di relazione tra spesa e
% umero di visite è virtualmente pari a 0 (8.5999e-170). Rifiutiamo
% decisamente l'ipotesi nulla di assenza di relazione tra numero di visite
% al supermercato e ammontare della spesa effettuata.


%% Calcolo manuale della matrice dei p-values
% 3) Calcolare manualmente la matrice dei p-values

% n= sample size
[n,p]=size(X);
% 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
Pval1=tcdf(abs(Testt),n-2,'upper')*2;
% Aggiungo i valori 1 sulla diagonale principale della matrice Pval1
Pvalchk=Pval1 +eye(p);
% Controllo che il massimo delle differenze in valore assoluto tra le due
% implementazioni sia pari ad un numero vicino a zero 
max(max(abs(Pval-Pvalchk)))

%% Esercizio aggiuntivo
% In un campione bivariato di 12 pezzi da un processo produttivo il
% coefficiente di correlazione è risultato pari a 0.54.
% Verificare che il p-value teorico del test (nell'ipotesi di ipotesi
% alternativa bilateale) è pari a 0.0699.
% Se invece di un solo campione di 12 pezzi noi disponessimo di 100000
% campioni di 12 pezzi, se è vera l'ipotesi nulla di assenza di
% correlazione tra le due variabili su quanti campioni ci aspettiamo di un
% ottenere un valore del coefficiente di correlazione superiore a 0.54 in
% modulo?
n=12;
r=0.54;
testR=sqrt(n-2)*r/sqrt(1-r^2);
pvalTeorico=tcdf(testR,n-2,'upper')*2;
% Se avessimo potuto disporre invece di un solo campione di 12 pezzi di
% 100000 campioni di 12 pezzi la frequenza relativa di campioni in cui il
% coefficiente di correlazione lineare sarebbe risultato superiore di 0.54
% in modulo sarebbe stata del 7 per cento circa.

%% Verifica empirica della distribuzione T di student con n-2 gdl
% 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 con quelli teorici. Utilizzare 100000 simulazioni ed un
% valore di n a piacere. 
% 
% Utilizzare i quantili che seguono.
% 
% quant=[0.01 0.05 0.10:0.1:0.9 0.95 0.99 0.999];
% 
n=12;
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.01 0.05 0.10:0.1:0.9 0.95 0.99 0.999]';
% Calcolo dei quantili empirici
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')

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

%% Verifico empiricamente il numero di volte in cui |r|>0.54
% Dato che il p-value teorico è pari al 7 per cento circa ci attendiamo che
% che su 100000 repliche nel 7 per cento di casi (7000) il valore del
% coefficiente di correlazione sia superiore a 0.54
pvalEmpirico=sum(abs(Rall)>r)/nsimul;
disp('p-value teorico calcolato tramite la distribuzione T(n-2)')
disp(pvalTeorico)
disp('p-value teorico calcolato estraendo 100000 campioni')
disp(pvalEmpirico)



%% Valore soglia per la significatività di r al variare di n
% 7) Rappresentare graficamente la soglia oltre il quale il valore di rxy è
% ritenuto significativo al livello del 5% in funzione della numerosità
% campionaria
% 
figure
nn=5:1000;
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')