%% LABORATORIO 5 - METODI STATISTICI PER LA BIOINGEGNERIA
%% Canale B - A.A. 2024/2025
%% Docente: Enrico Longato

%% Parte 1 (già fatta in un laboratorio precedente)
clc        
clear all  
close all  

%% Caricamento dati
load IVGTT_SS

%% Pulizia dati: individuare NaN, valori infiniti, valori negativi
% Qui come maschera; indifferente se preferite usare gli indici col find
i_mask_nan = isnan(data);
i_mask_inf = isinf(data);
i_mask_neg = data < 0;

%% Stampare il numero di dati di ciascuna categoria nell'intera tabella
% Calcolo con sum (potevo fare length del vettore di indici col find)
n_nan_tot = sum(i_mask_nan(:)); % oppure sum(sum(i_mask_nan)
n_inf_tot = sum(i_mask_inf(:));
n_neg_tot = sum(i_mask_neg(:));

disp(' ') % Riga vuota per cosmesi
disp('Analisi su tutta la tabella')
disp(['  - Numero di valori mancanti: ' num2str(n_nan_tot)])
disp(['  - Numero di valori infiniti: ' num2str(n_inf_tot)])
disp(['  - Numero di valori negativi: ' num2str(n_neg_tot)])

%% Pulizia dati: sostituire valori infiniti e valori negativi con NaN
disp(' ')
disp('Poiché non ci sono valori infiniti, né negativi, non faccio nulla.')
% Risposta valida all'esame.

%% Parte 2 (perfezionata da un laboratorio precedente: v. "!NEW!")
% Uso un ciclo for "per esercizio"; all'esame vanno bene soluzioni anche 
% meno (o più) raffinate
figure
for i=1:length(labels) % Scandisco entrambe le variabili

    %% Calcolo di statistiche descrittive per ogni variabile
    
    % Seleziono la colonna
    variable_curr = data(:, i);

    % Individuo e conto i dati mancanti
    i_nan_curr = isnan(variable_curr);
    n_nan_curr = sum(i_nan_curr);

    % Tolgo i dati mancanti (altrimenti la risposta è NaN, v. lab 3)
    variable_curr = variable_curr(~i_nan_curr);

    % Calcolo media, sd, mediana e IQR
    mean_curr = mean(variable_curr);
    std_curr = std(variable_curr);
    median_curr = median(variable_curr);
    iqr_curr = iqr(variable_curr);

    % Skeweness e curtosi
    skewness_curr = skewness(variable_curr);
    kurtosis_curr = kurtosis(variable_curr);

    % !NEW! Test di gaussianità (Lilliefors perché è quello che c'è in
    % MATLAB)
    [h_gauss, p_gauss, kstat_gauss, critval_gauss] = ... % ... = a capo riga
                               lillietest(variable_curr);

    % "Per abitudine" faccio il display sulla command window tutto insieme.
    % Si poteva fare man mano
    disp(' ')
    disp(['La variabile ' labels{i} ' (' units{i} ') ha:'])
    disp(['  - Numero di valori mancanti: ' num2str(n_nan_curr)])
    disp(['  - Media (senza NaN): ' num2str(mean_curr)])
    disp(['  - SD (senza NaN): ' num2str(std_curr)])
    disp(['  - Mediana (senza NaN): ' num2str(median_curr)])
    disp(['  - IQR (senza NaN): ' num2str(iqr_curr)])
    disp(['  - Skeweness (senza NaN): ' num2str(skewness_curr)])
    disp(['  - Curtosi (senza NaN): ' num2str(kurtosis_curr)])
    disp(['  - p-value del Lilliefors test: ' num2str(p_gauss)])

    %% Rappresentazione grafica delle distribuzioni
    % Boxplot nella prima colonna della figura
    subplot(2, 3, (i-1)*3 + 1) 
    boxplot(variable_curr)
    xlabel(labels{i})
    ylabel(units{i})
    
    % Istogrammi nella seconda colonna 
    subplot(2, 3, (i-1)*3 + 2) 
    N_bins = 20; % Per l'istogramma
    [counts_curr, centres_curr] = hist(variable_curr, N_bins);
    bar(centres_curr, counts_curr)
    xlabel(units{i})
    ylabel('Frequenze assolute (count)')
    subtitle(labels{i})

    % !NEW! QQ plot nella terza colonna
    subplot(2, 3, (i-1)*3 + 3)
    qqplot(variable_curr)
    subtitle(labels{i})
    % N.B.: Non complichiamo troppo le cose: le ordinate avrebbero le
    % unità di misura, le ascisse no (domanda per casa: perché?);
    % all'esame va bene se il QQ-plot ha le etichette di default (tranne
    % eventualmente il titolo)
    
end

% !NEW!
% Lascio il disp() dei miei commenti (sono del tenore delle risposte
% d'esame; come si vede, non serve sperticarsi in grandi dissertazioni)
disp(' ')
disp('Sulla base dei plot e delle statistiche descrittive stampate sopra, possiamo dire che')
disp('  - Nessuno dei due istogrammi sembra, a prima vista, incompatibile con la gaussianità (no bimodalità, no evidenti asimmetrie)')
disp('  - Nessuno dei due test ha rifiutato l''ipotesi nulla => Purtroppo, non possiamo dire nulla.')
disp('  - I due QQ-plot giacciono sostanzialmente su una retta; le code sono di scarsissimo rilievo: le due variabili sembrano ragionamente distribuite, ciascuna, come una gaussiana')
disp('  - Skewness e curtosi sono, in ambo i casi, vicine a 0 e a 3, rispettivamente. Anche questo è compatibile con la presunta gaussianità delle due varibaili.')
disp('      => In conclusione, nonostante il test statistico non abbia dato un risultato utile, gli altri indicatori sono tali da permetterci di assumere che le due variabili siano gaussiane.')

%% Parte 3
%% t-test a due campioni a due code
glucose_elderly = data(~isnan(data(:, 1)), 1);
glucose_young = data(~isnan(data(:, 2)), 2);

% Stessa sintassi di ttest, ma con un argomento in più (la seconda
% variabile) e un 2 alla fine del nome della funzione
[h, p, ci, stats] = ...
    ttest2(glucose_elderly, glucose_young, 'Vartype', 'unequal');

% Mostro i risultati e faccio una deduzione
disp(' ') % cosmesi
disp('L''ipotesi di gaussianità era stata verificata ed assunta per buona')
disp('H0 era che le due medie fossero uguali')
disp('Alfa era stato lasciato di default a 0.05.')
disp('Il t-test a due campioni indipendenti era a due code.')
disp(['Il p-value è risultato ' num2str(p)])
disp(['Non posso rifiutare H0, quindi non posso dedurre nulla sulla' ...
      ' diversità o meno tra le due medie di popolazione.'])