## Copper-nickel alloy, page 291 in Givens & Hoeting book

y <- c(127.6, 124.0, 110.8, 103.9, 101.5, 130.1, 
       122.0, 92.3, 113.1, 83.7, 128.0, 91.4, 86.2)
x <- c(0.01, 0.48, 0.71, 0.95, 1.19, 0.01, 0.48, 
       1.44, 0.71, 1.96, 0.01, 1.44, 1.96)

ndata <- length(y)

mylm <- function(y, x){

    fit <- lm(y ~ x)
    coefs <- coef(fit)
    theta <- as.numeric(coefs[2]/coefs[1])
    return(theta)
}

originalFit <- mylm(y,x)

B <- 10000

# decide which observation to include 
#(at random and with replacement)
obs.idx <- matrix(sample(1:ndata, ndata*B, replace=T), 
                  ncol=B, byrow=F)
# apply the regression model to each data subset
bootest <- sapply(1:B, function(i){
        mylm(y[obs.idx[,i]], x[obs.idx[,i]])
    }
)

## equivalent to using a for-loop
#bootest2 = rep(0, B)
#for(i in 1:B){
#  bootest2[i] =  mylm(y[obs.idx[,i]], x[obs.idx[,i]])
#}


library(MASS)
truehist(bootest, prob=TRUE, 
    ylab="Density", xlab=expression(theta), 
    col="lightblue")
abline(v=originalFit, col=2, lwd=3)

# determine bias
mean(bootest - originalFit)

# bias corrected estimate
originalFit - mean(bootest - originalFit)


# compute 95% CI based on percentile method.
quantile(bootest, prob=c(0.025, 0.975))