## jh 2016.09.04

## estimate parameters in Celsius to Fahrenheit conversion
## using 4 perfect (or imperfectly measured) C values

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

set.seed(  1234567 ) # personalize, 
                   ### so each student gets difefrent results 

n.estimates <- 1000

C <- seq(14,20,2) ; C

F <- 32 + 1.8*C ; cbind(C,F)

lm(F ~ C)$coefficients



par(mfrow=c(2,5),mar = c(5,5,1,1)) # 2 panels in graphics
 
 
### MEASUREMENT ERRORS ONLY IN F (i.e. on Y axis) ###########

sigmaF <- 2
sigmaC <- 0

beta.hat.if.error.in.F <- rep(NA,n.estimates)

## ************   fill in the blank 2 lines below   

for (dataset in 1:n.estimates) {

 f <- F + rnorm(4,mean = 0,sd = ****** ) ## errors in F
 beta.hat.if.error.in.F[dataset] <- lm(f ~ C)$coefficients[2]
 if (dataset <= 4) {
 	plot(C,f, pch=19,
           ylim=range(F)+c(-3,3),
           xlim=range(C)+c(-3,3),cex=1.5,
           main=paste(toString(dataset),
                      " of ", toString(n.estimates),sep=""))
    lines(C,lm(f ~ C)$fitted.values)
    points(C,F, pch=19,col="green",cex=0.75)
    lines(C,F, col="green",lwd=0.5)
 	
 } 
 
}

hist(beta.hat.if.error.in.F,breaks=30,
  xlim=c(0.5,4), main=paste("All ", toString(n.estimates),sep=""))

summary(beta.hat.if.error.in.F)
c( mean(beta.hat.if.error.in.F),
   sd(beta.hat.if.error.in.F), 
   var(beta.hat.if.error.in.F) )

se.mean.beta.hat.if.error.in.F <- sd(beta.hat.if.error.in.F)/
                      sqrt(n.estimates)

mean(beta.hat.if.error.in.F) + 
 c(-1.96, 1.96) * se.mean.beta.hat.if.error.in.F 



### MEASUREMENT ERRORS IN C (i.e. on X axis) ###########

sigmaF = 0; 
sigmaC=1;

beta.hat.if.error.in.C <- rep(NA,n.estimates)

## ************   fill in the blank 2 lines below

for (dataset in 1:n.estimates) {

 c <- C + rnorm(4,mean = 0,sd = ****** ) ## errors in C

 beta.hat.if.error.in.C[dataset] <- lm(F ~ c)$coefficients[2]
 if (dataset <= 4) {
 	 plot(c,F, pch=19,
           ylim=range(F)+c(-3,3),
           xlim=range(C)+c(-3,3),cex=1.5,
           main=paste(toString(dataset),
                      " of ", toString(n.estimates),sep=""))
     points(C,F, pch=19,col="green",cex=0.75)
     lines(C,F, col="green",lwd=0.5)
     lines(c,lm(F ~ c)$fitted.values)
 } 
 
 
}

hist(beta.hat.if.error.in.C,breaks=75,
 xlim=c(0.5,4), 
  main=paste("All ", toString(n.estimates),sep=""))

summary(beta.hat.if.error.in.C)

c( mean(beta.hat.if.error.in.C), 
   sd(beta.hat.if.error.in.C), 
   var(beta.hat.if.error.in.C) )

se.mean.beta.hat.if.error.in.C <- sd(beta.hat.if.error.in.C)/
sqrt(n.estimates)

mean(beta.hat.if.error.in.C) + 
 c(-1.96, 1.96) * se.mean.beta.hat.if.error.in.C