# Central Limit theorem
 # Written for Advanced Statistics, Econ. NTU
 # Jin-Lung Lin, April 2007
set.seed(12345678);
nobs=500;
nrep=500;
n0=10;
p=.25
muS=n0*p;
muU=0.5;
std0=sqrt(n0*p*(1-p));
rS=numeric(nrep);
rU=numeric(nrep);
rU1=numeric(nrep);
rU2=numeric(nrep);
rU3=numeric(nrep);
for (iobs in 20:nobs)
{
  for ( irep in 1:nrep)
  {
    S=rbinom(iobs,n0,p);
    U=runif(iobs);
    rS[irep]=sqrt(iobs)*(mean(S)-muS)/sqrt(var(S));
    rU[irep]=sqrt(iobs)*(mean(U)-muU)/sqrt(var(U));
   } 
 op=par(mfcol=c(2,2))
 plot(rS); 
 hist(rS,ylim=c(0,1),xlim=c(-5,5),breaks=15,prob=T,col="magenta");
 curve(dnorm(x),-5,5,add=TRUE);
 mtext(paste("# of observations=",iobs),side=3);
 plot(rU); 
 hist(rU,ylim=c(0,1),xlim=c(-5,5),breaks=15,prob=T,col="magenta");
 curve(dnorm(x),-5,5,add=TRUE);
 mtext(paste("# of observations=",iobs),side=3);
 par(op)
}

for (iobs in 20:nobs)
{
  for ( irep in 1:nrep)
  {
     U=runif(iobs);
     rU[irep]=sqrt(iobs)*(mean(U)-muU)/sqrt(var(U));
     rU1[irep]=mean(U);
     rU2[irep]=iobs*(mean(U)-muS);
     rU3[irep]=sqrt(iobs)*mean(U)/sqrt(var(U));    
   } 
 op=par(mfcol=c(2,2))
 hist(rU,ylim=c(0,1),xlim=c(-5,5),breaks=15,prob=T,col="magenta");
 mtext(paste("Mean of", "# of observations=",iobs),side=3);
 hist(rU1,ylim=c(0,1),xlim=c(-5,5),breaks=15,prob=T,col="magenta");
  mtext(paste("mean(X)", "# of observations=",iobs),side=3);
  hist(rU2,breaks=15,prob=T,col="magenta");
 mtext(paste("T * (mean(U)-mu)", "# of observations=",iobs),side=3);
  hist(rU3,breaks=15,prob=T,col="magenta");
 mtext(paste("T * (mean(U))", "# of observations=",iobs),side=3);
 par(op)
}





for ( i in 10:nobs)
 { x=rnorm(i);
   hist(x,ylim=c(0,1),xlim=c(-5,5),prob=T);
   curve(dnorm(x),-5,5,add=TRUE);
   mtext(paste("# of observations=",i),side=3);
   }