blue1 = function(n){
y = rep(0,1000)
for(i in 1:1000) y[i] = sum(runif(n)<.0588)/n
hist(y,nclass=30,main=paste("n = ",as.character(n)))
}

par(mfrow=c(3,3))
for(n in c(1,5,10,25,50,75,100,200,500)) blue1(n)

x = (runif(10000)<.0588)
y = cumsum(x)/c(1:10000)
par(mfrow=c(1,1))
plot(c(1,10000),c(0,.1),xlab="n",ylab="x-bar",type="n",main="plot of x-bar")
lines(c(1:10000),y,col=rainbow(10)[1])
lines(c(1:10000),.0588+2*0.23525/sqrt(1:10000),lty=2,col=rainbow(10)[4])
lines(c(1:10000),.0588-2*0.23525/sqrt(1:10000),lty=2,col=rainbow(10)[4])
abline(h=.0588,lty=2,col=rainbow(10)[7])