clear 
close all
clc
%% A) Visualizzazione dei dati

% Carico i dati
load('data_IVGTT_SS.mat')
%%
% Estraggo le variabili in due vettori colonna
Gss_elderly = data(:,1);
Gss_young = data(:,2);

% Eliminare eventuali valori negativi, NaN o inf perché non fisiologici

if(any(isnan(Gss_elderly)) || any(isinf(Gss_elderly)) || any(Gss_elderly<0))
    disp('Ci sono NaN, Inf o valori negativi: non sono fisiologici e vanno eliminati!') 
    Gss_elderly(isnan(Gss_elderly)) = [];
    Gss_elderly(isinf(Gss_elderly)) = [];
    Gss_elderly(Gss_elderly<0) = [];
end

if(any(isnan(Gss_young)) || any(isinf(Gss_young)) || any(Gss_young<0))
    disp('Ci sono NaN, Inf o valori negativi: non sono fisiologici e vanno eliminati!')  
    Gss_young(isnan(Gss_young)) = [];
    Gss_young(isinf(Gss_young)) = [];
    Gss_young(Gss_young<0) = [];
end


% Calcolare media, standard deviation, mediana e moda delle distribuzioni delle due variabili e visualizzare i valori nella command window 
mean_Gss_elderly = mean(Gss_elderly);
mean_Gss_young = mean(Gss_young);

sigma_Gss_elderly = std(Gss_elderly);
sigma_Gss_young = std(Gss_young);

median_Gss_elderly =  median(Gss_elderly);
median_Gss_young = median(Gss_young);

mode_Gss_elderly = mode(Gss_elderly);
mode_Gss_young = mode(Gss_young);

disp('----------------------------------------------') 
disp(['Media di Gss_elderly = ', num2str(mean_Gss_elderly)]);
disp(['Media di Gss_young = ', num2str(mean_Gss_young)]);

disp(['STD di Gss_elderly = ', num2str(sigma_Gss_elderly)]);
disp(['STD di Gss_young = ', num2str(sigma_Gss_young)]);

disp(['Mediana di Gss_elderly = ', num2str(median_Gss_elderly)]);
disp(['Mediana di Gss_young = ', num2str(median_Gss_young)]);

disp(['Moda di Gss_elderly = ', num2str(mode_Gss_elderly)]);
disp(['Moda di Gss_young = ', num2str(mode_Gss_young)]);


%% B) Verifica l'assunzione di gaussianità del test di Student:
 
% Controllare la distribuzione delle due variabili con i boxplot

figure(1)
boxplot(data,labels)

% Controllare la distribuzione delle due variabili con gli istogrammi

figure(2)
subplot(1,2,1)
histogram(Gss_elderly,'Normalization','pdf');
title('Pdf G-elderly')
ylabel('Frequenze normalizzate')
xlabel('Gss')
hold on

subplot(1,2,2)
histogram(Gss_young,'Normalization','pdf')
title('Pdf G-young')
ylabel('Frequenze normalizzate')
xlabel('Gss')
hold on

[    ,   ] = ksdensity(    );
pdf_elderly = normpdf(    ,    ,    );

subplot(1,2,1)
% Plot kernel-density
plot(   ,   ,'-g','LineWidth',2)
% Plot pdf from normpdf (gaussian distribution)
plot(   ,   ,'-b','LineWidth',2)
hold off

% ks-density + pdf gaussiana + istogramma
[   ,   ] = ksdensity(   );
pdf_young = normpdf(   ,    ,    );

subplot(1,2,2)
% Plot kernel-density
plot( , ,'-g','LineWidth',2)
% Plot pdf from normpdf (gaussian distribution)
plot( , ,'-b','LineWidth',2)
hold off

% Calcolare kurtosis e skewness e visualizzare i valori nella command window
kurtosis_Gss_young   =  kurtosis(Gss_young);
kurtosis_Gss_elderly =  kurtosis(Gss_elderly);
disp('----------------------------------------------') 
disp(['Kurtosis di Gss_young = ', num2str(kurtosis_Gss_young)]);
disp(['Kurtosis di Gss_elderly = ', num2str(kurtosis_Gss_elderly)]);

skewness_Gss_young   =  skewness(Gss_young);
skewness_Gss_elderly =  skewness(Gss_elderly);
disp(['Skewness di Gss_young = ', num2str(skewness_Gss_young)]);
disp(['Skewness di Gss_elderly= ', num2str(skewness_Gss_elderly)]);

% Verifica la gaussianita di Gss_young e Gss_elderly con lillie-fors test
[h,p]=lillietest(Gss_young);
[h,p]=lillietest(Gss_elderly);


%% C) Calcolare la statistica t (utilizzare le formule nelle slides)

% Gruppo 1 - young 
% Gruppo 2 - elderly

N1=numel(Gss_young);
N2=numel(Gss_elderly);

mu1=mean_Gss_young;
mu2=mean_Gss_elderly;

sigma1=sigma_Gss_young;
sigma2=sigma_Gss_elderly;

sp2 = 


t_obs =


%% D) Calcolare la distribuzione di Student nel range x = [-10:0.01:10] (tpdf)

x = [-10:.01:10];
df = N1+N2-2; %gradi di libertà
t_tstudent = tpdf(        );
figure(3)
plot(x,t_tstudent)
title('Distribuzione di student')
xlabel('t')

%% E) Calcolare il pvalue e verificare se l'ipotesi nulla (H0) è accettata o rifiutata considerando un livello di significatività α = 5%. 

% Si calcoli il valore dell'integrale a partire da t
area_coda_dx =           ;
pvalue = 

disp('-----------------------')
disp(['Pvalue = ',num2str(pvalue)])

%% F) Calcolare il valore critico (TH) corrispondente ad alpha significativo testando i valori del range x a partire dal valore 1  


%% G) Confrontando la t_obs con il valore critico (TH), verificare se lGipotesi nulla (H0) è accettata. (BONUS: visualizzare i valori di tobs e TH nel plot della distribuzione di Student)

% Plot
figure(4)


%% H) Confrontare i risultati calcolati a mano con i parametri h e p restituiti dalla function matlab ttest2.

[h,p,ci,stats] = ttest2(                 );

disp('-----------------------')
disp(['Pvalue restituito da ttest2 = ',num2str(p)]);
disp(['Risultato del test restituito da ttest2 = ',num2str(h)]);
disp(['Statistica t restituita da ttest2 = ',num2str(stats.tstat)]);