%% Testo dell'esercizio
% 1) Caricare i dati contenuti nel file Firm.xlsx in formato table dentro
% una variabile denominata X utilizzando i dati presenti nella colonna A
% come rownames della table.
% 
% 2) Calcolare il massimo il minimo e la mediana di ogni variabile
% (quantitativa)
% 
% 3) Trasformare le variabili peso e altezza rispettivamente in kg e cm.
% Creare una nuova table esattamente uguale alla precedente ma con le
% variabili peso e altezza in Kg e cm. 
% Le formule di trasformazione sono:
% Hcm= (H inches) *2.54
% Wkg= (Weight pounds)/2.2046;
% Inserire come nomi delle nuove variabili Wkg e Hcm.
% Risolvere il quesito utilizzando i metodi che seguono:
% 3A Transform into cm and kg using function table2array. Put the output in
% a table called YA
% 3B Transform into cm and kg extracting the contents of the variables
% using the dot notation and manipulating directly. Put the output in
% a table called YB
% 3C Transform into cm and kg extracting the variables using curly brackets
% notation and the corresponding number of the column. Put the output in
% a table called YC
% 3D Transform into cm and kg extracting the variables using curly brackets
% and the corresponding names of the variables. Put the output in
% a table called YD.
% Put the results using the different method in a separate cell and make
% sure that YA = YB = YC = YD (hint: use function isequal)
%
% Put the variables height in cm, weight in kg and wage in a matrix of
% doubles called HWW (of size nx3), where n is the number of rows of the
% datasets. For the questions below use matrix HWW
%
% 4) Calcolare medie, deviazioni standard, varianze, minimo e massimo delle 
% variabili peso, altezza e salario
% 
% 5) Calcolare la matrice nx3 degli scostamenti e degli scostamenti standardizzati 
% del peso, dell'altezza e del salario
% 5A) utilizzando l'espansione implicita tramite il vettore riga che
% contiene le medie ed il vettore riga che contiene le standard deviation
% 5B) utilizzando un loop (ciclo for)
% 5C) utilizzando la funzione di MATLAB zscore
% 
% 6) Rappresentare graficamente gli scostamenti standardizzati delle righe che 
% iini:ifin (inserendo come coordinate sull'asse X i
% nomi delle righe). Assegnare alla variabile iini il valore 1 ed alla
% variabile ifin il valore 10
% 
% 7) Calcolare gli indici di asimmetria e di curtosi (v. i lucidi e/o
% l'help delle funzioni MATLAB skewness e kurtosis per le formule di calcolo).
% 7A) Utilizzare l'implementazione manuale (inserendo direttamente in MATLAB la
% formula di calcolo)
% 7B) Confrontare il risultato con quello che si ottiene facendo
% riferimento direttamente alle funzioni MATLAB skewness e kurtosis
%
% 8) Calcolare (e interpretare) i quantili che seguono [0.025 0.25 0.50 0.75 
% 0.975]  della variabile retribuzione (wage)
% 
% 9) Rappresentare graficamente (e commentare) il boxplot della variabile 
% altezza (per tutte le unitā e poi distinto per sesso)
% 
% 10) Calcolare la distribuzione di frequenza della variabile sesso (gender)  
% e fornirne una rappresentazione grafica adeguata
% 
% 11) Analizzare la variabile retribuzione (wage).
% Assegnare la table iniziale ad una nuova variabile denominata Ycor 
% (sempre in formato table). Costruire il boxplot della variabile retribuzione 
% facendo riferimento ad Ycor. Modificare l'unitā che presente una retribuzione 
% superiore a 10000 con il valore 1916.26. Costruire il boxplot della variabile 
% retribuzione modificata.
% 
% 12)  Costruzione del grafico ad istogrammi (utilizzando la funzione histogram) 
% della retribuzione modificata, prima utilizzando il numero di classi di default 
% poi un numero di classi pari a 15 e poi utilizando le classi
% 1500-2000, 2000-2500, ..., 3500-4000.
% Tramite la funzione subplot inserire i tre grafici di cui sopra in 3
% pannelli separati
%
% 
% 13) Trovare le frequenze associate alle classi 1500-2000, 2000-2500, ...,
% 3500-400 chiamando la funzione histogram tramite l'argomento di output
% ossia h = histogram(___)
% 
% 14) Calcolare le frequenze tramite un ciclo for
% Calcolare le frequenze utilizzando le classi
% '<1500' '1500-2000' '2000-2500'   '2500-3000' '3000-3500' '3500-4000' '>4000'
% in maniera alternativa utilizzando un ciclo for combinato con l'istruzione if.
% Inserire le frequenze trovate dentro una table che contenga come nomi
% delle righe le etichette di cui sopra.
%
% 15) Costruzione della tabella di contingenza tra le variabili Education  
% e Gender utilizzando la funzione crosstab. Mostrare l'output in formato table.
% 
% 16) Una volta terminato script salvare il file in formato .m ed in
% formato .mlx con il nome analisi_preliminari_out
% 
% 
% 
% *SOLUZIONE*
% 
% 
%% Read table inside MATLAB
% 1) Caricare i dati contenuti nel file Firm.xlsx in formato table dentro
% una variabile denominata X utilizzando i dati presenti nella colonna A
% come rownames della table.

X=readtable('Firm.xlsx','ReadRowNames',true);

% Se i dati sono caricati senza l'opzione 'ReadRowNames',true per
% aggiungere i nomi delle righe occorre operarare come segue
%{
    % read withtout rownames
    Xchk=readtable('Firm.xlsx');
    % Modo 1 
    % Insert the content of column 1 (transformed into cell) as rownames
    Xchk.Properties.RowNames=table2cell(Xchk(:,1));
    % Modo2 
    % Extract the contents of column 1 and put it in the RowNames property
    % of the table
    Xchk.Properties.RowNames=Xchk{:,1};
    % Delete first column (because now it is rownames)
    Xchk=Xchk(:,2:end);
%}


% Per un tutorial sui vari tipi di dati presenti in MATLAB si veda
% https://www.mathworks.com/help/matlab/data-types_data-types.html

%% 2) Calcolare il massimo il minimo e la mediana di ogni variabile (quantitativa)
summary(X)

%% 3) Trasformare le variabili peso e altezza rispettivamente in kg e cm.
% Creare una nuova table esattamente uguale alla precedente ma con le
% variabili peso e altezza in Kg e cm. 
% Le formule di trasformazione sono:
% Hcm= (H inches) *2.54
% Wkg= (Weight pounds)/2.2046;
% Inserire come nomi delle nuove variabili Wkg e Hcm.

%% 3A Transform into cm and kg using function table2array
%  Transform variable height into cm
Hcm=table2array(X(:,7))*2.54;

% Transform variable weight into kg
Wkg=table2array(X(:,8))/2.2046;

% Create an n-by-2 array which contains  
% height in cm and the weight in kg 
HW=[Hcm Wkg];

% Create an n-by-2 table which contains  
% height in cm and the weight in kg 
HWTable=array2table(HW,'VariableNames',{'Heightcm', 'Weightkg'});

% YA = table di dimensione nx8 che contiene nelle colonne 7 e 8 l'altezza
% in cm ed il peso in kg
YA=[X(:,1:6) HWTable];

%% 3B Transform into cm and kg extracting the contents of the variables using the dot notation and manipulating directly 
% Di seguito si illustra un modo alternativo per trasformare i dati in una
% table senza passare attraverso le funzioni table2array e array2table
YB=X;
YB.Height=YB.Height*2.54;
YB.Weight=YB.Weight/2.2046;
YB.Properties.VariableNames(7)={'Heightcm'};
YB.Properties.VariableNames(8)={'Weightkg'};
% Controllo che la table YA sia esattamente uguale alla table YB
disp('Check equality of YA and YB')
isequal(YA,YB)

%% 3C Transform into cm and kg extracting the variables using curly brackets and the corresponding number of the column
YC=X;
% I dati dalla table vengono estratti direttamente in formato double
% tramite le parentesi graffe
YC{:,7}=YC{:,7}*2.54;
YC{:,8}=YC{:,8}/2.2046;
YC.Properties.VariableNames(7)={'Heightcm'};
YC.Properties.VariableNames(8)={'Weightkg'};
disp('Check equality of YA and YC')
isequal(YA,YC)

%% 3D Transform into cm and kg extracting the variables using curly brackets and the corresponding names of the variables
YD=X;
% I dati dalla table vengono estratti direttamente in formato double
% tramite le parentesi graffe
% Le colonne della table vengono estratte direttamente passando il loro
% nome
YD{:,'Height'}=YD{:,'Height'}*2.54;
YD{:,'Weight'}=YD{:,'Weight'}/2.2046;
YD.Properties.VariableNames('Height')={'Heightcm'};
YD.Properties.VariableNames('Weight')={'Weightkg'};
disp('Check equality of YA and YD')
isequal(YA,YD)

%% Create matrix HWW  which contains H W and X.Wage
% Put the variables height in cm, weight in kg and wage in a matrix of
% doubles called HWW (of size nx3), where n is the number of rows of the
% datasets. For the questions below use matrix HWW
HWW= [HW X.Wage];

%% 4. Calcolo delle medie, deviazioni standard, varianze, minimo e massimo

% Calcolo delle medie
AvHWW=mean(HWW);

% Calcolo dei sigma
% Osservazione: std ha come denominatore n-1 (stimatore corretto di sigma)
stdHWW=std(HWW);

% Calcolo delle varianze
% Osservazione: var ha come denominatore n-1 (stimatore corretto della varianza)
varHWW=var(HWW);

% Calcolo dei mimini
minHW=min(HWW);

% Calcolo dei massimi
maxHWW=max(HWW);

%% 5 Calcolare la matrice nx3 degli scostamenti e degli scostamenti standardizzati 
% del peso, dell'altezza e del salario
% 5A) utilizzando l'espansione implicita tramite il vettore riga che
% contiene le medie ed il vettore riga che contiene le standard deviation
% 5B) utilizzando un loop (ciclo for)
% 5C) utilizzando la funzione di MATLAB zscore

%% Calcolo delle matrici degli scostamenti e degli scostamenti standardizzati del peso e dell'altezza
% 5A) Risoluzione tramite espansione implicita

% Matrice degli scostamenti dalla media utilizzando l'espansione implicita
% Compute deviations from means
% HWW_Tilde = Matrice di dimensione n-by-3 che contiene gli scostamenti
% dalla media
% Ogni colonna della matrice HeiWeiTilde presenta media 0
% HWW č d dimensione 107x3
% AvHWW č di dimensioni 1x3
% Tramite la cosiddetta espansione implicita MATLAB (a partire dalla
% versione 2017A) capisce che il vettore di dimensione 1x3 deve applicarsi
% ad ogni riga della matrice HWW 
HWW_Tilde=HWW-AvHWW;

% Matrice degli scostamenti standardizzati
% Compute zscore of height and weight
% Calcolo degli scostamenti standardizzati
% Ogni colonna della matrice HW_Z presenta media zero e varianza 1
% Si noti che il vettore stdHW presenta dimensione 1x3 la matrice HW_Tilde
% (matrice degli scostamenti dalla media) presenta dimensione nx3.
% La notazione A ./ B (right array division) divide ciascun elemento di A per 
% il corrispondente elemento di B. Tramite l'espansione implicita (dalla
% release 2017A) Matlab sa che deve dividere ogni riga (di dimensione
% 1x3) della matrice HW_Tilde per il corrispondente elemento del vettore
% stdHW (di dimensione 1x3)
HWW_Z=HWW_Tilde./stdHWW;


%% Zscores computed using a loop
% 5B) utilizzando un loop (ciclo for)
% Calcolo degli scostamenti standardizzati tramite un ciclo loop
% Inizializzazione della matrice che conterrā gli scostamenti
% standardizzati
n=size(HWW,1);

% Viene inizializzata la matrice di dimensione nx3 che conterrā gli
% scostamenti standardizzati calcolati tramite ciclo for
HWW_ZWithLoop=zeros(n,3);
for j=1:3
   HWW_ZWithLoop(:,j)=(HWW(:,j)-mean(HWW(:,j)))/std(HWW(:,j));
end
% Controllo che la matrice degli scostamenti standardizzati ottenuta sia la stessa che č stata ottenuta con
% i metodi precedenti
isequal(HWW_ZWithLoop,HWW_Z)

%% Compute zscores using MATLAB function zscore
% 5C) utilizzando la funzione di MATLAB zscore
% Calcolo gli scostamenti standardizzati per il peso e l'altezza tramite la
% funzione zscore
HW_ZCHK=zscore(HWW);

% Controlliamo che gli scostamenti standardizzati calcolati con i due
% metodi siano gli stessi
isequal(HWW_Z,HW_ZCHK)

%% 6. Rappresentazione grafica degli scostamenti standardizzati

% Rappresentazione tramite grafico a barra degli scostamenti standardizzati
% delle unitā statistiche iini:ifin (inserendo come coordinate sull'asse X
% i nomi delle righe)

% Osservazione: nell'help della funzione bar si legge che 
% "One way to indicate categories for your bars is to specify X as a
% categorical array." E' necessario quindi trasformare i nomi delle righe
% (che costituiscono un tipo di dati cell) in un tipo di dati categorico
% tramite la funzione categorical
% Osservazione:  
% class(X.Properties.RowNames(1:10)) riporta
%    'cell'

iini=1;
ifin=10;
bar(categorical(X.Properties.RowNames(iini:ifin)),HWW_Z(iini:ifin,:))

% Aggiunta della legenda al grafico
legend({'Altezza' 'Peso'})
title('Scostamenti standardizzati')
%% 7. Indici di asimmetria (skewness)
% 7) Calcolare gli indici di asimmetria e di curtosi (v. i lucidi e/o
% l'help delle funzioni MATLAB skewness e kurtosis per le formule di calcolo).
% 7A) Utilizzare l'implementazione manuale (inserendo direttamente in MATLAB la
% formula di calcolo)
% 7B) Confrontare il risultato con quello che si ottiene facendo
% riferimento direttamente alle funzioni MATLAB skewness e kurtosis
%

% Osservazione: X.^2 eleva al quadrato ogni singolo elemento della matrixe
% X

%% Indice di asimmetria
% Implementazione manuale facendo riferimento direttamente alle formule

% s2 = vettore riga che contiene la media dei quadrati degli scostamenti
% dalla media (varianze non corrette)
s2 = mean(HWW_Tilde.^2); 
% Osservazione: s2 = var(HW_Tilde,1) =var(HW,1)

% m3= vettore riga che contiene la media delle potenze terze degli
% scostamenti dalla media
m3 = mean(HWW_Tilde.^3);
% AS = indice di asimmetria non corretto
AS = m3 ./ s2.^(1.5);
% AScorr = stimatore corretto della asimmetria .
AScorr = AS .* sqrt((n-1)./n) .* n./(n-2);

% Controllo che l'indice di asimmetria che č stato implementato č esattamente
% uguale a quello ottenuto tramite la funzione di MATLAB skewness
AScorrchk=skewness(HWW,0);
disp('Controllo uguaglianza tra le due implementazioni dell''indice di asimmetria')
disp(isequal(AScorr,AScorrchk))
   
   
%% Indici di curtosi
% Implementazione manuale facendo riferimento direttamente alle formule
% Indice di curtosi non corretto
% m4= vettore riga che contiene la media delle potenze quarte degli
% scostamenti dalla media
m4 = mean(HWW_Tilde.^4); 
k = m4 ./ s2.^2;

% kcor= indice di curtosi corretto
kcor = ((n+1)*k - 3*(n-1)) .* (n-1)/((n-2)*(n-3)) + 3;

% Controllo che l'indice di curtosi che č stato implementato č esattamente
% uguale a quello ottenuto tramite la funzione di MATLAB kurtosis
kcorchk=kurtosis(HWW,0);
disp('Controllo uguaglianza tra le due implementazioni dell''indice di curtosi')
disp(isequal(kcor,kcorchk))

%% QUANTILI
% 8) Calcolare (e interpretare) i quantili che seguono [0.025 0.25 0.50 0.75 
% 0.975]  della variabile retribuzione (wage)
% 
% Calcolo dei quantili richiesti per la variabile retribuzione (Wage)
quant=[0.025 0.25 0.50 0.75 0.975];
qq=quantile(X.Wage,quant);

% Istruzione alternative passando attraverso table2array
% qq=quantile(table2array(X(:,5)),quant);
% Istruzione alternativa senza passare attraverso table2array
% qq= quantile(X{:,5},quant)

%% BOXPLOTS
% 9) Rappresentare graficamente (e commentare) il boxplot della variabile 
% altezza (per tutte le unitā e poi distinto per sesso)
% 

Y=YA;
% Calcolo boxplot per la variabile altezza
boxplot(YA.Heightcm)

ylabel('Altezza in cm')


%% Boxplot per la variabile altezza distinto per sesso
% Istruzione alternativa senza passare attraverso table2array

% Istruzione alternativa chiamando le colonne richieste con il loro nome nella table
boxplot(Y.Heightcm,Y.Gender)
% Istruzioni alternative
% boxplot(Y{:,7},Y{:,3})
% boxplot(table2array(Y(:,7)),table2array(Y(:,3)))

xlabel('Sesso')
ylabel('Altezza in cm')

% Modo alternativo di estrarre la variabile altezza
% boxplot(table2array(Y(:,'Heightcm')))
% boxplot(table2array(Y(:,'Heightcm')),table2array(Y(:,'Gender')))


%% Analisi della variabile Sesso (Gender)
% 10) Calcolare la distribuzione di frequenza della variabile sesso
% (gender) e fornirne una rappresentazione grafica adeguata
% 

% Distribuzione di frequenza della variabile Gender
out=tabulate(Y.Gender);

% Si noti che out č un oggetto di tipo cell

% Costruzione del pie chart (grafico a torta)
% Osservazione: la funzione pie necessita di un input numerico 
pie([out{:,2}])
% Metodo alternative per estrarre gli elementi numerici nella colonna 2
% della table
% pie(cell2mat(out(:,2)))

% Vengono aggiunte come labels le etichette della colonna 1
pie([out{:,2}],out(:,1))


% Grafico a torta della variabile Gender con label personalizzate
label={'Uomo' 'Donna'};
pie([out{:,2}],label)

% Metodo alternativo per costruire il grafico a torta trasformando il cell
% array Y.Gender in una variabile categorica e passarandola direttamente
% alla funzione pie
pie(categorical(Y.Gender))

%% Analisi della variabile retribuzione
% 11) Analizzare la variabile retribuzione (wage).
% Assegnare la table iniziale ad una nuova variabile denominata Ycor 
% (sempre in formato table). Costruire il boxplot della variabile retribuzione 
% facendo riferimento ad Ycor. Modificare l'unitā che presente una retribuzione 
% superiore a 10000 con il valore 1916.26. Costruire il boxplot della variabile 
% retribuzione modificata.

% Il boxplot di questa variabile rivela un valore anomalo molto grande L'unitā 
% 4 č un chiaro valore anomalo!
boxplot(Y.Wage)

% Assegnare la table iniziale ad una nuova variabile denominata Ycor 
% (sempre in formato table)
Ycor=Y;
% Modificare l'unitā che presente una retribuzione 
% superiore a 10000 con il valore 1916.26.
Ycor.Wage(Ycor.Wage>10000) = 1916.26;
% Osservazione: l'istruzione Ycor.Wage>10000 produce un vettore Booleno che
% presente un true (1) in corrispondenza delle righe per cui la condizione
% č vera e un faso (0) quando la condizione č falsa.
% In questo esempio l'istruzione č vera solo in corrispondenza della riga
% 4.

%% Costruzione dei grafici ad istogrammi della retribuzione
% 12)  Costruzione del grafico ad istogrammi (utilizzando la funzione histogram) 
% della retribuzione modificata, prima utilizzando il numero di classi di default 
% poi un numero di classi pari a 15 e poi utilizando le classi
% 1500-2000, 2000-2500, ..., 3500-4000.
% Tramite la funzione subplot inserire i tre grafici di cui sopra in 3
% pannelli separati
% L'istruzione figure serve per aprire una nuova finestra grafica
figure

subplot(3,1,1)
histogram(Ycor.Wage)
title('Grafico a istogrammi utilizzando il numero di classi di default')

subplot(3,1,2)
% Costruzione del grafico ad istogrammi specificando 15 classi
histogram(Ycor.Wage,15)
title('Grafico a istogrammi utilizzando 15 classi')

subplot(3,1,3)
% Costruzione del grafico ad istogrammi specificando gli intervalli
% richiesti per le classi
classes=1500:500:4000;
histogram(Ycor.Wage,classes)
title('Grafico a istogrammi utilizzando le classi di ampiezza 500 partendo da 1500')


%% Recupero delle frequenze visualizzate nel grafico 
% 13) Trovare le frequenze associate alle classi 1500-2000, 2000-2500, ...,
% 3500-400 chiamando la funzione histogram tramite l'argomento di output
% ossia h = histogram(___)

% trovare le frequenze associate alle classi di cui sopra
figure
h=histogram(Ycor.Wage,1500:500:4000);
disp('Frequenze associate alle diverse classi')
disp(h.Values)
% Aggiunta delle etichette
% text((1750:500:3750)',(h.Values+3)',num2str(h.Values'))

%% 14) Calcolare le frequenze tramite un ciclo for
% Calcolare le frequenze utilizzando le classi
% '<1500' '1500-2000' '2000-2500'   '2500-3000' '3000-3500' '3500-4000' '>4000'
% in maniera alternativa utilizzando un ciclo for combinato con l'istruzione if.
% Inserire le frequenze trovate dentro una table che contenga come nomi
% delle righe le etichette di cui sopra.

% Estraggo la variabile Wage
Wage=Ycor.Wage;

% distfreq inizializzo il vettore che conterrā le 7 frequenze in
% corrispondenza delle 7 classi richieste
distfreq=zeros(7,1);
for i=1:length(Wage)
    if Wage(i)<1500
        distfreq(1)=distfreq(1)+1;
    elseif Wage(i)>=1500 && Wage(i) <2000
        distfreq(2)=distfreq(2)+1;
    elseif Wage(i)>=2000 && Wage(i) <2500
        distfreq(3)=distfreq(3)+1;
    elseif Wage(i)>=2500 && Wage(i) <3000
        distfreq(4)=distfreq(4)+1;
    elseif Wage(i)>=3000 && Wage(i) <3500
        distfreq(5)=distfreq(5)+1;
    elseif Wage(i)>=3500 && Wage(i) <4000
        distfreq(6)=distfreq(6)+1;
    else
        distfreq(7)=distfreq(7)+1;
    end
end

% Creo una cell array of strings che contiene i nomi delle classi
Nomiclassi={'<1500' '1500-2000' '2000-2500' ...
    '2500-3000' '3000-3500' '3500-4000' '>4000'};

% l'array distrfreq viene trasformato in table
% Le opzioni 'VariableNames' e 'RowNames' consentono di personalizzare i
% nomi delle righe e delle colonne
Freqtable=array2table(distfreq,'VariableNames',{'Frequenze'},'RowNames',Nomiclassi);
disp(Freqtable)

%% Costruzione tabelle di contingenza
% 15) Costruzione della tabella di contingenza tra le variabili Education  
% e Gender utilizzando la funzione crosstab. Mostrare l'output in formato table.
 [conting,chi2,p,labels] =crosstab(Y.Education,Y.Gender);
 [r,c]=size(conting);
% Costruzione della tabella di contingenza in formato table (con i nomi
% delle righe e delle colonne)
 ContTable=array2table(conting,'RowNames',labels(1:r,1),'VariableNames',labels(1:c,2));

disp(ContTable)