mu=10
sigma=2
n=300
X<- rnorm(n,mu,sigma) # A sample of size n from N(10,1^2).
hist(X,freq=FALSE,xlim=c(mu-3*sigma,mu+3*sigma),ylim=c(0,1.2*dnorm(mu,mu,sigma)))
lines(seq(mu-3*sigma,mu+3*sigma,length=100),dnorm(seq(mu-3*sigma,mu+3*sigma,length=100),mu,sigma),col="blue")
X
barX=mean(X)
cat("Based on this sample, sample mean barX is", barX, "\n",
    "where the true population mean is", mu)

# Let us simulate the sample distribution of the sample mean
# I will simulate M values of the sample mean based on
# M different random samples, each of which is of size n
n=30
M=1000
SampleMean=c()
for(m in 1:M)
{
  X<- rnorm(n,mu,sigma)
  SampleMean=c(SampleMean,mean(X))
}
hist(SampleMean,freq=FALSE,xlim=c(mu-3*sigma,mu+3*sigma),ylim=c(0,4),ylab="Sampling distribution of the sample mean", xlab="")
legend("topright",legend=paste("n=",n))
lines(seq(mu-3*sigma,mu+3*sigma,length=100),dnorm(seq(mu-3*sigma,mu+3*sigma,length=100),mu,sigma),col="blue")
segments(mu,0,mu,4,col="red")

par(mfrow=c(2,2))
par(mar=c(3,5,4,.1))
for(n in c(4, 16,25,100))
{
  M=1000
  SampleMean=c()
  for(m in 1:M)
  {
    X<- rnorm(n,mu,sigma)
    SampleMean=c(SampleMean,mean(X))
  }
  hist(SampleMean,freq=FALSE,xlim=c(mu-3*sigma,mu+3*sigma),ylim=c(0,2),ylab="Sampling distribution of the sample mean", xlab="")
  legend("topright",legend=paste("n=",n))
  lines(seq(mu-3*sigma,mu+3*sigma,length=100),dnorm(seq(mu-3*sigma,mu+3*sigma,length=100),mu,sigma),col="blue")
  segments(mu,0,mu,4,col="red")
}


# Sample variance

par(mfrow=c(2,2))
par(mar=c(3,5,4,.1))
for(n in c(4, 16,25,100))
{
  M=1000
  SampleVariance=c()
  for(m in 1:M)
  {
    X<- rnorm(n,mu,sigma)
    SampleVariance=c(SampleVariance,var(X))
  }
  hist(SampleVariance,freq=FALSE,xlim=c(0,5),ylim=c(0,3),ylab="Sampling distribution of the sample variance", xlab="")
  legend("topright",legend=paste("n=",n))
  segments(sigma^2,0,sigma^2,3,col="red")
}