% Enrico Longato
% 9929778
% Canale B

% Il file deve iniziare come scritto sopra e avere un nome file come da
% consegna del testo d'esame.

% Pulizia workspace
clc
clear all
close all

% Caricamento dati
load simulazione_20250113

%% Esercizio 1
% Individuare e riportare il numero di valori mancanti e valori non
% fisiologici
N_nan = sum(isnan(data(:))); % Numero di NaN in tutta la tabella
N_neg = sum(data(:) < 0); % Numero di valori negativi in tutta la tabella

% Eliminare le colonne con >20% dati mancanti o non fisiologici
i_nan_or_neg_mask = isnan(data) | (data < 0); % Maschera dei valori NaN o negativi

perc_nan_by_col = 100 * sum(i_nan_or_neg_mask) / size(data, 1); % Percentuale di NaN/negativi per colonna

data_reduced = data(:, perc_nan_by_col<=20); % Tengo tutte le righe, le colonne con <= 20% NaN/negativi

% Eliminare le righe in cui siano rimasti valori NaN o non fisiologici
i_no_nans = ~any(isnan(data_reduced) | (data_reduced < 0), 2); % 2 = esamino ogni riga
data_reduced_filtered = data_reduced(i_no_nans, :); % Tengo le righe senza NaN e tutte le colonne della matrice ridotta

% Riportare le dimensioni
[dim_1, dim_2] = size(data_reduced_filtered);

%% Esercizio 2
% Verificare la gaussianità della prima colonna di data
age = data(:, 1);

figure % L'ordine dei pannelli è indifferente; l'importante sono le etichette e le unità di misura
subplot(121)
histogram(age)
title('Istogramma delle frequenze assolute')
xlabel('Età (anni)')
ylabel('Numero di osservazioni')

subplot(122)
qqplot(age)
title('QQ-plot')

age_skewness = skewness(age); % Da confrontare con il valore 0
age_kurtosis = kurtosis(age); % Da confrontare con il valore 3

[h_gauss, p_gauss, kstat_gauss, critval_gauss] =  lillietest(age);

% Test statistico per confrontare la media con il valore 46.6
target = 46.6;
[h, p, ci, stats] = ttest(age, target);

%% Esercizio 3
% Individuare X e Y
X = [data(:, [1, 5, 6]) ones(size(data, 1), 1)]; % Intercetta inclusa
Y = data(:, 7);

% Stimare "a mano" i parametri
beta_hat = (X'*X)\X'*Y;

% Intervallo di confidenza intorno al parametro relativo a "medication"
sigma2_hat = (sum((Y - X*beta_hat).^2))/(size(data, 1) - length(beta_hat)); % Varianza a posteriori

XtX_inv = inv(X'*X);
standard_errors = sqrt(diag(sigma2_hat*XtX_inv));

i_medication = 3; % Bisogna tenerne traccia "a mano"; dipende da dove era stata messa l'intercetta
beta_medication = beta_hat(i_medication); 
standard_error_medication = standard_errors(i_medication);

confidence_interval_medication = [beta_medication - 1.96 * standard_error_medication; ...
                                  beta_medication + 1.96 * standard_error_medication]; % andava bene anche usare +/-2

%% Esercizio 4
% Disegnare le tracce
figure
plot(signal_time, signal_data)
title('Segnali misurati')
xlabel('Tempo (minuti)')
ylabel('Concentrazione (mg/dL)')

% Utilizzare due volte l'algoritmo k-means con le specifiche date
N_replicates = 10; % Ripetizioni del K-means per non finire in ottimi locali
N_iter = 300; % Numero massimo di iterazioni dell'algoritmo per ogni replica

rng(42) % seed
[cluster_2, centroid_2] = kmeans(signal_data, 2, ...
                                 'Distance', 'sqeuclidean', ...
                                 'Replicates', N_replicates, ...
                                 'MaxIter', N_iter);
rng(42) % seed
[cluster_3, centroid_3] = kmeans(signal_data, 3, ...
                                 'Distance', 'sqeuclidean', ...
                                 'Replicates', N_replicates, ...
                                 'MaxIter', N_iter);

% Valori di silhouette per K = 2 o 3
silhouette_2 = mean(silhouette(signal_data, cluster_2));
silhouette_3 = mean(silhouette(signal_data, cluster_3));