par(mfrow=c(2,1)) 		# 2 by 1 array of plots
x <- seq(1,16,0.1)
y <- 0.6*dnorm(x,5,1)+0.4*dnorm(x,10,2)
plot(x,y,type="l")
x <- matrix(rep(0,60),6) 	# 6 by 10 matrix full of zeroes
n <- rep(0,6)
ss <- rep(0,60)
sm <- rep(0,60)
index <- 0			# 6 by 1 vector of zeroes
for(i in 1:6){
 n[i] <- 10^(i-1) 		# sample sizes are n=1, 10, 100, 1000, 10000, 100000
 for(j in 1:10){
  index <- index + 1
  sample.mean <- 0
  for(k in 1:n[i]){
   B <- rbinom(1,1,0.6) 	# 1 draw from bin(1,0.6)
   C1 <- rnorm(1,5,1)
   C2 <- rnorm(1,10,2)
   sample.mean <- sample.mean + B*C1+(1-B)*C2
  }
  x[i,j] <- sample.mean/n[i]
  ss[index] <- i
  sm[index] <- x[i,j]
 }
}
plot(ss,sm,main="10 sample means from sizes n=1,10,100,1000,10000,100000")