nreps=5000;
n=5;
mu=20;
sigma=1;
for i=1:nreps;
    tx=normrnd(mu,sigma,n,1);
    mvec(i)=mean(tx);
    svec(i)=var(tx);
    tvec(i)=(mvec(i)-mu)/(sqrt(svec(i)/n));
end

lowert=min(tvec);
uppert=max(tvec);
step=(uppert-lowert)/25;
xbr=lowert:step:uppert;
subplot(2,2,1);
hist(tvec,xbr);
hold on;
x=lowert+(uppert-lowert)*(0:1000)/1000;
yz=normpdf(x,0,1);
plot(x,nreps*step*yz);
title('T Statistic: Sampling distribution')
legend('N(0,1) pdf')
axis([lowert uppert 0 2500]);
hold off

subplot(2,2,2);
hist(tvec,xbr);
hold on;
yt=tpdf(x,n-1);
plot(x,nreps*step*yt);
title('T Statistic: Sampling distribution')
legend('Student(n-1) pdf')
axis([lowert uppert 0 2500]);
hold off

subplot(2,2,3);
hold on
plot(x,yt)
plot(x,yz,'green')
legend('N(0,1) pdf','Student(n-1)',2)
axis([-5 5 0 0.75]);
hold off


    


