%% Testo
% 1) Caricare  le variabili  "retta" e "affiliazione" presenti nel file di Excel UnivUSAanova.xlsx tramite la funzione xlsread. Denominare la variabile retta X e la variabile affiliazione group
% 2) rimuovere in maniera listwise le righe del dataset  che presentano valori mancanti o inammissibili in almeno una delle due variabili  (valori mancanti o inammissibili = valori che presentano stringa uguale a NaN nella variabile affiliazione e dati mancanti nella variabile retta).
% 3) Costruire la distribuzione di frequenze della variabile classificatoria group.
% 4) Analizzare la distribuzione della retta tramite boxplot nei 3 strati della variabile classificatoria ("Clero" "Stato"  "Privati").
% 5) Tramite la funzione di MATLAB anova1 testare l’ipotesi che le rette medie per iscriversi a scuole statali, private o del clero differiscano tra loro significativamente. Calcolare e commentare il p-value del test. Commentare il risultato del test.
% 6) Calcolare la devianza tra i gruppi (DEVtra) e nei gruppi (DEVnei). Implementare manualmente il test sull'uguaglianza delle medie e calcolare il relativo p-value
% 7) Calcolare la media campionaria e l'intervallo di confidenza (utilizzando il livello di confidenza specificato in una variabile denominata conflev) della media dell'universo della variabile X (retta pagata) per ciascuno dei 3 strati. Calcolare gli estremi dell'intervallo di confidenza utilizzando sia la formula per grandi campioni, sia la formula per piccoli campioni
%
%
%
% 8) Rappresentare graficamente i 3 intervalli di confidenza tramite la funzione errorbar
% 9) [QUESITO AVANZATO] Costruire una funzione (denominata ConfInt) con due argomenti di input obbligatori  e due  facoltativi. Il primo argomento di input obbligatorio deve essere un vettore di lunghezza n di dati numerici. Il secondo argomento di input deve essere una cell o un vettore numerico di lunghezza n contenente i valori di una variabile qualitativa classificatoria. Il primo argomento facoltativo deve essere un booleano che, se è true, in output mi deve apparire su grafico la rappresentazione grafica degli intervalli di confidenza della media dei dati numerici per ogni modalità (strato) della variabile classificatoria. Il valore di default del primo argomento facoltativo di input deve essere false. Il secondo argomento di input facoltativo deve essere un scalare di tipo double compreso nell'intervallo (0 1) che specifica l'ampiezza dell'intervallo di confidenza delle medie dell'universo nei diversi strati. Se l'utente non specifica il secondo argomento di input facoltativo utilizzare come default un livello di confidenza uguale a 0.95.
% Osservazione: testare che l'utente fornisca in input gli argomenti richiesti obbligatori con le specifiche richieste e fornire un messaggio di errore se questo non avviene.
% La funzione deve avere tre argomenti di output.
% Il primo argomento di output deve essere un array di dimensione k x 3 (dove k è il numero di livelli delle variabile classificatoria) contenente nella prima colonna gli estremi inferiori degli intervalli di confidenza, nella seconda colonna le medie campionarie e nelle terza colonna gli estremi superiori degli intervalli di confidenza.
% Il secondo argomento deve essere una table di dimensione kx3 contenente le informazioni di cui sopra e i nomi appropriati delle righe e delle colonne
%  Il terzo argomento di output deve essere una table di dimensione 1 x 2 contenente la quota di varianza (devianza totale) spiegata da DEVnei e la quota di devianza spiegata da DEVtra.
% 10) Utilizzando la variabile retta X e la variabile affiliazione group chiamare la funzione Confint creata al punto precedente, mostrando il grafico in output e impostando a 99.99% il valore dell'intervallo di confidenze delle medie dell'universo nei tre strati.
% 11) Calcolare delle medie, varianze e intervalli di confidenza per ogni gruppo tramite la funzione grpstats.
% 12) Rappresentare graficamente gli intervalli di confidenza al 99% delle medie dei strati tramite la chiamata alla funzione grpstats
%% Soluzione
%% 1) Caricamento dati da file di Excel con xlsread

[X, group, raw] = xlsread('UnivUSAanova.xlsx','UnivUSA','K2:L251');

%
%% 2) Rimozione valori mancanti e inammissibili

NaNindexesGroup=strcmp(group,'NaN');
% conta i valori di group che sono uguali alla stringa NaN
% sum(NaNindexesGroup)

[a,indexmissX]=rmmissing(X);
% conta i missing di X
% sum(indexmissX)

% NaNindexes = vettore Booleano che presenta true se una riga è missing in
% X oppure se presenta la stringa NaN dentro group
NaNindexes=NaNindexesGroup | indexmissX;
% conta i missing totali
sum(NaNindexes)
% Opera sui dati non mancanti
X=X(~NaNindexes);
group=group(~NaNindexes);

%
%% 3) Distribuzione di frequenze della variabile classificatoria


tabulate(group)

%% 4) Analisi preliminare tramite boxplot

boxplot(X,group);

%% 5) Test ANOVA di uguaglianza delle medie


% Analisi ANOVA
[p,anovatab,stats]=anova1(X,group)

%
% Commento: il p-value è molto basso di conseguenza rifiuto decisamente l'ipotesi di uguaglianza delle tre medie nell'universo.
%
%% 6) e 7) Calcolo manuale della devianza nei gruppi e della devianza tra i gruppi e degli intervalli di confidenza

% unigroup=cell che contiene le etichette dei diversi gruppi 
unigroup=unique(group);

% inizializzazione di devNEI e devFRA  
devNEI=0;
devFRA=0;
n=length(X);
ngroups=length(unigroup);
% conflev = scalare associato all'ampiezza dell'intervallo di confidenza
conflev=0.99;
% zalpha= quantile 0.995, ossia valore che lascia alla sua destra il 0.5
% per cento dei valori in una v.c. normale standardizzata
% quant= 0.995 
quant=1-(1-conflev)/2;
zalpha=norminv(quant);
Mglobal=mean(X);

% Inizializzazione della matrice ConfN
% Prima colonna di Conf = estremo inferiore intervallo di confidenza del
% gruppo j (utilizzando la formula per grandi campioni)
% Seconda colonna di Conf = Media del gruppo j 
% Terza colonna di Conf = estremo superiore intervallo di confidenza del
% gruppo j (utilizzando la formula per grandi campioni)
ConfN=zeros(ngroups,3);

% La matrice ConfT presenta le stesse dimensioni e lo stesso significato
% della matrice confN ma utilizza la formula dell'intervallo di confidenza
% per piccoli campioni.
ConfT=ConfN;

for j=1:ngroups
    booj=strcmp(group,unigroup(j));
    % nj = numerosità del gruppo j
    nj=sum(booj);
    % Xj= vettore di lunghezza nj che contiene le unità appartenenti al
    % gruppo j
    Xj=X(booj);
    % varj = varianza (corretta) delle unità appartenenti al gruppo j
    varj=var(Xj);
    % stdj= standard deviation (s corretto) delle unità appartenenti ai
    % gruppo j
    stdj=sqrt(varj);
    
    devNEI=devNEI+varj*(nj-1);
    % Mj=media del gruppo j
    Mj=mean(Xj);
    devFRA=devFRA+(Mj-Mglobal)^2 *nj;
    
    % Costruzione dell'intervallo di confidenza per le unità appartenenti
    % al gruppo j (utilizzando la formula per grandi campioni)
    
    

%

    ConfN(j,:)=[Mj-zalpha*stdj/sqrt(nj) Mj  Mj+zalpha*stdj/sqrt(nj)];
    % Costruzione dell'intervallo di confidenza per le unità appartenenti
    % al gruppo j (utilizzando la formula per piccoli campioni)
    

%

    talpha=tinv(quant,nj-1);
    ConfT(j,:)=[Mj-talpha*stdj/sqrt(nj) Mj  Mj+talpha*stdj/sqrt(nj)];
    
end
varFRA=devFRA/n;
varNEI=devNEI/n;

% test ANOVA implementazione manuale
numANOVA=devFRA/(ngroups-1);
denANOVA=devNEI/(n-ngroups);
testANOVA=numANOVA/denANOVA;
pvaltestANOVA=fcdf(testANOVA,ngroups-1,n-ngroups,'upper');

%% 8) Rappresentazione grafica intervalli di confidenza


% Rappresentazione degli intervalli di confidenza tramite la funzione
% errorbar
errorbar((1:ngroups)',ConfN(:,2),(ConfN(:,3)-ConfN(:,2)),'o')
% L'istruzione che segue aggiunge nella parte alta del grafico le etichette definite dentro unigroup
text([0.2; 0.5; 0.8],1.05*ones(ngroups,1),unigroup,'Units','normalized','HorizontalAlignment','Center')
% definisce il limite minimo e massimo dell'asse x
xlim([0.5 3.5])
set(gca,'XTickLabel','')
ylabel(['Intervalli di confidenza al ' num2str(100*conflev) '% della retta da pagare'])

%% 10) Utilizzando la variabile retta X e la variabile affiliazione group chiamare la funzione Confint creata al punto precedente, mostrando il grafico in output e impostando a 99.99% il valore dell'intervallo di confidenze delle medie dell'universo nei tre strati.
%

[ConfN,Conftable,devtable]=Confint(X,group,true,0.9999)


%% Calcolo delle medie, varianze e intervalli di confidenza per ogni gruppo direttamente tramite la funzione grpstats

% Utilizzando le due variabili X e group creo una table di dimensione nx2
% con nomi delle variabili retta e affliazione
Xtable=table(X,group,'VariableNames',{'retta','affiliazione'});
X1=grpstats(Xtable,'affiliazione',{'mean','std','var','meanci'});
% Mostre le medie le std le varianze e gli intervalli di confidenza al 95
% per cento basati utilizzando la formula per piccoli campioni
disp(X1)

% Osservazione: se grpstats viene chiamato con 'Alpha',valore di alpha è possibile controllare l'ampiezza dell'intervallo di confidenza

X1=grpstats(Xtable,'affiliazione',{'mean','std','var','meanci'},'Alpha',0.01);
% Mostre le medie le std le varianze e gli intervalli di confidenza al 99
% per cento basati utilizzando la formula per piccoli campioni
disp(X1)



%
%% Rappresentazione grafica intervalli di confidenza tramite gprstats

%% Rappresentazione grafica intervalli di confidenza al 99 per cento
% Dall'help di grpstats
% grpstats(X,group,alpha) plots the means of the groups of data in the
% vector or matrix X determined by the values of the grouping variable,
% group. The grouping variable values are on the horizontal plot axis. Each
% group mean has 100×(1 – alpha) per cent confidence intervals.
% Osservazione: gli intervalli di confidenza utilizzano i quantili della
% v.v. T di Student (e non quelli della normale)
grpstats(X,group,0.01);