%% 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)