nreps=5000;
n=10;
mu=20;
sigma=1;
for i=1:nreps;
    mvec(i)=mean(normrnd(mu,sigma,n,1));
end
lower=min(mvec);
upper=max(mvec);
step=(upper-lower)/25;
xbr=lower:step:upper;
subplot(1,2,1);
hist(mvec,xbr);
hold on;
x=lower+(upper-lower)*(0:1000)/1000;
y=normpdf(x,mu,sigma/sqrt(n));
plot(x,nreps*step*y);
title('Mean: Sampling distribution')
axis([lower upper 0 750]);
hold off

zvec=(mvec-mu)/(sigma/sqrt(n));
lower=min(zvec);
upper=max(zvec);
step=(upper-lower)/25;
xbr=lower:step:upper;
subplot(1,2,2);
hist(zvec,xbr);
hold on;
x=lower+(upper-lower)*(0:1000)/1000;
y=normpdf(x,0,1);
plot(x,nreps*step*y);
title('Z Statistic: Sampling distribution')
legend('N(0,1) pdf')
axis([lower upper 0 750]);
hold off



    


