clear all
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

[fks_elderly,xi_elderly] = ksdensity(Gss_elderly);
pdf_elderly = normpdf(xi_elderly,mean_Gss_elderly,sigma_Gss_elderly);

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

% ks-density + pdf gaussiana + istogramma
[fks_young,xi_young] = ksdensity(Gss_young);
pdf_young = normpdf(xi_young,mean_Gss_young,sigma_Gss_young);

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

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


skewness_Gss_elderly =  skewness(Gss_elderly);
skewness_Gss_young   =  skewness(Gss_young);

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


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



%% 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 = (((N1-1)*sigma1^2)+ ((N2-1)*sigma2^2))/ ...
                  (N1+N2-2);

t_obs =        (mu1-mu2)/ ...
       (sqrt(sp2)*sqrt(1/N1+1/N2));

%% 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(x,df);
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 = (1-tcdf(t_obs,df));

pvalue = 2*area_coda_dx;
disp('-----------------------')
disp(['Pvalue = ',num2str(pvalue)])


% alternativa con trapz
% mask_tobs = x>t_obs;
% area_coda_dx = trapz(x(mask_tobs),t_tstudent(mask_tobs));


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

idx = find(x==1);

for i = idx:length(x)
    temp_alpha_mezzi = trapz(x(i:end),t_tstudent(i:end));
    temp_alpha = 2*temp_alpha_mezzi;
    if (temp_alpha<0.05)
        TH = x(i);
        disp(['Soglia (Threshold) = ', num2str(TH)]);
        break
    end
end

%% 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)
plot(x,t_tstudent,'k')
hold on
plot(t_obs,0,'b*','LineWidth',10)
plot(TH,0,'g*','LineWidth',10)
plot([t_obs,t_obs],[0,0.4],'b-','LineWidth',2)
plot([TH,TH],[0,t_tstudent(i)],'g-','LineWidth',2)
xlim([-4,4])
title('Distribution di t','FontSize',18)
xlabel('t','FontSize',18)
legend("Distribuzione di Student", "Statistica t", "Valore critico TH")
hold off

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

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

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)]);