%Traditional versus Piaget 
clear all
close all
% randn('state',0)
% trad = round(normrnd(84.1, 9,[1 16]))
% piaget = round(normrnd(88.1, 15,[1 14]))

trad =[80    69    85    87    74    85    95    84  ...
         87    86    82    91    79   100    83    85]; 
  
piaget =[ 100    89    87    76    93    68    99  ...
                  100    78   99   100   74    76    97];

X1bar = mean(trad) %84.5

X2bar = mean(piaget)  %88.2857
s1 = std(trad)  %7.3847
s2 = std(piaget)  %11.6515
n1= 16; n2 = 14;

 %=========================
 F = s1^2/s2^2  %F =0.4017
 pval1 = 2*fcdf(F, n2-1,  n1 -1) %pval1 =0.1063
 % decide sigma_1^2 = sigma_2^2 test by
 % pooled standard deviation

 sp = sqrt(  ( (n1-1)*s1^2 + (n2-1)*s2^2)/(n1 + n2 - 2) )
 %sp =9.6044
 %
 t = (X1bar - X2bar)/(sp * sqrt(1/n1 + 1/n2)) %t=-1.0771
 %
 pval = tcdf(t, n1 + n2 - 2) %pval=0.1453
 %
 tcrit = tinv(0.05, n1 + n2 - 2)  %tcrit=-1.7011
 %===============================================================
 %==============================================================
 
 disp('Let''s NOT assume equal variances')
 vt = var(trad);
 vp = var(piaget);
 T = (X1bar-X2bar)/sqrt(var(trad)/n1+var(piaget)/n2)
 df = (vt/n1+vp/n2)^2/((vt/n1)^2/(n1-1)+(vp/n2)^2/(n2-1))
 pvalue = tcdf(T,df)