# each of you should use a different seed to start the random number sequence 

sigmaF <- 1
sigmaC <- 1

estimates <- 1000

C <- seq(14,20,2)

F <- 32 + 1.8*C

lm(F ~ C)$coefficients


beta.hat.c <- rep(NA,estimates)

# jh's slower way

for (dataset in 1:estimates) {

 f <- F + rnorm(4,mean = 0,sd = sigmaF)

 beta.hat.c[dataset] <- lm(f ~ C)$coefficients[2]
 
}

hist(beta.hat.c)

summary(beta.hat.c)

c( mean(beta.hat.c), sd(beta.hat.c), var(beta.hat.c) )

se.mean.beta.hat.c <- sd(beta.hat.c)/sqrt(estimates)

mean(beta.hat.c) + c(-1.96, 1.96) * se.mean.beta.hat.c