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

%% Parte 1 
clc        
clear all  
close all  

%% Carichiamo i dati
load ELSA % Osservare il workspace dopo il caricamento

%% Ispezionando visivamente i dati, verifichiamo se sembrino continui o 
%% discreti e se eventuali valori negativi potrebbero essere plausibili

% La consegna richiedeva ispezione visiva. Possiamo stabilire se una
% variabile assume valori continui anche in maniera programmatica,
% verificando che la distanza tra due valori consecutivi sia sempre intera
% NB: potrebbero esserci casi particolari in cui questo algoritmo non 
% funziona, ma sono preoccupiamocene

disp(' ') % Salto una riga per cosmesi
for i_col=1:size(elsa, 2)
    
    % Colonna/variabile corrente
    col_curr = elsa(:, i_col);

    % Riordino per avere i valori consecutivi adiacenti
    col_curr = sort(col_curr);

    % Calcolo il vettore delle differenze prime (cioè la differenza tra gli
    % elementi in posizione k+1 e k
    diff_curr = diff(col_curr);

    % Sottraggo alle differenze il loro valore approssimato all'intero più
    % vicino (se fossero interi, la differenza, ove sensata, sarebbe 0)
    diff_with_nearest_integer_curr = diff_curr - round(diff_curr);

    % Verifico se almeno un elemento è diverso da zero
    is_not_integer = any(diff_with_nearest_integer_curr ~= 0);

    % Rispondo
    if is_not_integer
        disp(['La variabile #' num2str(i_col) ' = ' elsa_labels{i_col} ...
              ' è continua'])
    else
        disp(['La variabile #' num2str(i_col) ' = ' elsa_labels{i_col} ...
              ' non sembra continua (verificare che non sia' ...
              ' semplicemente campionata con passo 1)'])
    end

end

% La plausibilità dei valori negativi è una constatazione fattibile solo "a
% buon senso" e/o sulla base della conoscenza di dominio.
% Nel caso specifico, nessuna delle variabili ha senso possa assumere
% valori negativi (che senso avrebbero un'età negativa o una pressione
% negativa?)

%% Valori mancanti (NaN) nell'intera tabella (naive)
% Prima possibilità usando una maschera booleana
i_nan_mask = isnan(elsa(:)); % Srotolo 
n_nan = sum(i_nan_mask); % Sommo gli 1 (logical true)

% Seconda possibilità usando find
i_nan = find(isnan(elsa)); % o elsa(:), tanto è sempre l'indice lineare
n_nan_2 = length(i_nan); % Lunghezza degli indici a cui c'era il nan

%% Valori registrati come infinito (Inf) nell'intera tabella (naive)
% Prima possibilità usando una maschera booleana
i_inf_mask = isinf(elsa(:)); % Srotolo 
n_inf = sum(i_inf_mask); % Sommo gli 1 (logical true)

% Seconda possibilità usando find
i_inf = find(isinf(elsa)); % o elsa(:), tanto è sempre l'indice lineare
n_inf_2 = length(i_inf); % Lunghezza degli indici a cui c'era l'inf

%% Valori negativi non fisiologici nell'intera tabella (naive)
% Prima possibilità usando una maschera booleana
i_negative_mask = elsa(:) < 0; % Srotolo 
n_negative = sum(i_negative_mask); % Sommo gli 1 (logical true)

% Seconda possibilità usando find
i_negative = find(elsa < 0); % o elsa(:), tanto è sempre l'indice lineare
n_negative_2 = length(i_negative); % Lunghezza degli indici a cui c'era l'inf

%% Stampiamo a schermo i risultati dell'analisi naive sull'intera matrice
disp(' ') % Salto una riga, per cosmesi
disp('Risultati dell''analisi naive sull''intera tabella:')
% '' è la escape sequence dell'apice
disp(['    - Numero di NaN = ' num2str(n_nan)])
disp(['    - Numero di Inf = ' num2str(n_inf)])
disp(['    - Numero di valori negativi = ' num2str(n_negative)])

% NB: L'analisi sull'intera tabella è poco interessante; queste
% considerazioni iniziano ad avere senso a livello di singola variabile e/o
% di singolo soggetto (da cui "naive")

%% Provo a calcolare medie e varianze
% (reminder: di default, sono la media e la varianza di ciascuna colonna)

mean_with_nans = mean(elsa);
var_with_nans = var(elsa);

% Solo il primo valore non è NaN. Scopriremo a breve che è perché c'è
% almeno un nan in ciascuna delle altre colonne!

%% Parte 2

%% Assegnare ai valori "inf" e negativi della matrice il valore "nan"
% Ragionamento: giacché sono valori implausibili, tanto vale considerarli 
% come mancanti

% Indici degli elementi che saranno nan (i nan di prima e il resto)
i_nan_new = isnan(elsa) | isinf(elsa) | (elsa < 0);

% Setto a nan (indicizzazione con maschera)
elsa(i_nan_new) = nan;

% NB alcuni lo erano già: poco male

%% Calcolare la percentuale di nan per ogni variabile
% Calcoliamo il numero di nan per colonna
nan_counts = sum(isnan(elsa));

% Trasformiamo in percentuale
nan_percentages = 100 * nan_counts / size(elsa, 1);

%% Rimuovere le colonne con >20% dati mancanti
elsa_reduced = elsa(:, nan_percentages <= 20);

%% Parte 3

%% Rimuovere i soggetti che hanno, in almeno una delle variabili rimaste, 
%% un valore mancante
% Sommo per righe con sum(..., 2) e guardo quando i nan sono più di 0
i_at_least_one_nan = sum(isnan(elsa_reduced), 2) > 0;
% Alternativamente, potevo fare any(isnan(elsa_reduced), 2)

% Filtro
elsa_reduced_filtered = elsa_reduced(~i_at_least_one_nan, :);

%% Applicare nuovamente l'operatore media/varianza
mean_without_nans = mean(elsa_reduced_filtered);
var_without_nans = var(elsa_reduced_filtered);
% Tutto come previsto

%% Parte 4

%% Copiare la matrice originale in una matrice elsa_imputed
elsa_imputed = elsa;


%% Sostituire ai valori mancanti il valore mediano dei dati non mancanti
%% della colonna a cui appartengono
for i_col = 1:size(elsa_imputed, 2)

    % Colonna corrente
    col_curr = elsa_imputed(:, i_col);
    
    % Indici dei nan
    i_nan_curr = isnan(col_curr);
    
    % Mediana dei non mancanti
    median_curr = median(col_curr(~i_nan_curr)); % Oppure nanmedian
    
    % Inserisco al posto giusto
    elsa_imputed(i_nan_curr, i_col) = median_curr;

end

%% Parte 5

%% Confrontare numerosità, media, varianza, mediana e istogrammi 
%% dei dati originariamente non mancanti di ciascuna colonna 
%% e di quelli dopo l'imputazione con la mediana

for i_col = 1:size(elsa_imputed, 2)
    
    % Colonna corrente nelle due matrici
    col_curr_original = elsa(:, i_col);
    col_curr_imputed = elsa_imputed(:, i_col);

    % Tolgo i nan dalla colonna originale
    col_curr_original = col_curr_original(~isnan(col_curr_original));

    % Numerosità
    N_original = length(col_curr_original);
    N_imputed = length(col_curr_imputed);

    % Media, varianza, mediana
    mean_original = mean(col_curr_original);
    median_original = median(col_curr_original);
    var_original = var(col_curr_original);

    mean_imputed = mean(col_curr_imputed);
    median_imputed = median(col_curr_imputed);
    var_imputed = var(col_curr_imputed);

    % Display
    disp(' ')
    disp(['Variabile ' elsa_labels{i_col}])
    disp(['Media originale vs. imputata: ' ...
           num2str(mean_original) ' vs. ' num2str(mean_imputed)])
    disp(['Mediana originale vs. imputata: ' ...
           num2str(median_original) ' vs. ' num2str(median_imputed)])
    disp(['Varianza originale vs. imputata: ' ...
           num2str(var_original) ' vs. ' num2str(var_imputed)])
    
    % Istogrammi
    N_bins = 20;
    [freq_original, centres_original] = hist(col_curr_original, N_bins);
    [freq_imputed, centres_imputed] = hist(col_curr_imputed, N_bins);

    % Plot
    figure

    subplot(121)
    bar(centres_original, freq_original)
    subtitle('Originale')
    ylabel('Frequenze assolute (count)')
    xlabel(elsa_units{i_col})

    subplot(122)
    bar(centres_imputed, freq_imputed)
    subtitle('Imputata')
    ylabel('Frequenze assolute (count)')
    xlabel(elsa_units{i_col})

    sgtitle(elsa_labels{i_col}) % Titolo sopra ai subplot

end

% NB: si poteva fare con un ciclo for su cell array di matrici

% NB 2: l'asse y cambia di scala (sono le frequenze assolute!) 
% NB 3: variabili che avevano tanti dati mancanti hanno un bin molto più
%       alto di prima dopo l'imputazione: è quello con dentro la mediana.