nreps=5000;
nx=5;
ny=9
mux=20;
muy=20;
sigmax=2;
sigmay=2;
for i=1:nreps;
    x=normrnd(mux,sigmax,nx,1);
    y=normrnd(muy,sigmay,ny,1);
    sigmapsq=((nx-1)*var(x)+(ny-1)*var(y))/(nx+ny-2);
    tvec(i)=(mean(x)-mean(y))/sqrt(sigmapsq*(1/nx+1/ny));
end

critvalue=tinv(0.975,nx+ny-2);

lowert=min(tvec);
uppert=max(tvec);
step=(uppert-lowert)/25;
xbr=lowert:step:uppert;
hist(tvec,xbr);
hold on;
x=lowert+(uppert-lowert)*(0:1000)/1000;
yt=tpdf(x,nx+ny-2);
plot(x,nreps*step*yt);
title('Two-sample T Statistic: Sampling distribution')
axis([lowert uppert 0 1500]);
legend('St(nx+ny-2) pdf')
line([critvalue critvalue],[0,1500],'Color','r');
line([-critvalue -critvalue],[0,1500],'Color','r');
legend('St(nx+ny-2) pdf','Critical Value')
t_rejectH0=length(tvec(abs(tvec)>critvalue))/nreps;
text(-2,1400,'Obs. Rejection Rate ')
text(-2,1300,num2str(t_rejectH0))

hold off



    


