n<-100
set.seed(101)
x<-runif(n,0,100)
be0<-2    #True value of beta0
be1<-1.5  #True value of beta1
ytrue<-2.0+1.5*x
plot(x,y,type='n')
points(x,ytrue,pch=19)
points(x,y,pch=19,col='red')

####################################
nreps<-5000

sigma<-20
beta.estimates<-matrix(0,nrow=nreps,ncol=2)
for(irep in 1:nreps){

	epsilon<-rnorm(n,0,sigma)
	y<-ytrue+epsilon
	beta.estimates[irep,]<-coef(lm(y~x))

}

hist(beta.estimates[,1],xlab=expression(hat(beta)[0]),main='',col='gray');box()
abline(v=be0,col='red',lty=2,lwd=2)
title(expression(paste('Histogram of 5000 estimates for ',beta[0])))

hist(beta.estimates[,2],xlab=expression(hat(beta)[1]),main='',col='gray');box()
abline(v=be1,col='red',lty=2,lwd=2)
title(expression(paste('Histogram of 5000 estimates for ',beta[1])))

plot(beta.estimates,pch=19,xlab=expression(hat(beta)[0]),ylab=expression(hat(beta)[1]))
abline(v=be0,h=be1,col='red',lty=2,lwd=2)
title(expression(paste('Scatterplot of 5000 pairs of estimates')))

####################################

X<-cbind(1,x)
XtXinv<-solve(t(X)%*%X)

Estimator.variance<-sigma^2*XtXinv

hist(beta.estimates[,1],xlab=expression(hat(beta)[0]),main='',col='gray',breaks=seq(-20,20,by=1));box()
abline(v=be0,col='red',lty=2,lwd=2)
title(expression(paste('Histogram of 5000 estimates for ',beta[0])))

xv<-seq(-20,20,by=0.01)
yv<-dnorm(xv,be0,sqrt(Estimator.variance[1,1]))  #The theoretical sampling distribution for beta[0]
lines(xv,yv*1*nreps,col='blue',lwd=2)


hist(beta.estimates[,2],xlab=expression(hat(beta)[1]),main='',col='gray',breaks=seq(1,2,by=0.025));box()
abline(v=be1,col='red',lty=2,lwd=2)
title(expression(paste('Histogram of 5000 estimates for ',beta[1])))

xv<-seq(1,2,by=0.001)
yv<-dnorm(xv,be1,sqrt(Estimator.variance[2,2]))  #The theoretical sampling distribution for beta[1]
lines(xv,yv*0.025*nreps,col='blue',lwd=2)