### R commands for sections 9.3
### G. Givens & J. Hoeting, Computational Statistics, Wiley, 2005.


## Example 9.8
temp <- read.table("alloy.dat", header=T)
data.size <- dim(temp)[1]
x <- rep(0, data.size)
y <- rep(0, data.size)
for(i in 1:data.size) {x[i] <- temp[i,1]; y[i] <- temp[i,2];}
# use linear regression to fit the model Y = beta0 + beta1*x + epsilon
temp <- lm(y~x)             # linear regression of y on x
summary(temp)
beta0 <- temp$coef[1]       # estimate for beta0
beta1 <- temp$coef[2]       # estimate for beta1
theta <- beta1/beta0        # estimate of beta1/beta1
theta

B0 <- B1 <- 300
# bootstrapping the cases
set.seed(2)
n <- length(x)              # size of original sample
R0star <- rep(0, B0)        # R0(X^*,F^hat) based on bootstrapped samples
R1star <- rep(0, B0)        # R1(X^*, F^hat) based on double bootstrapped samples
tempt <- proc.time()        # record system time
for(i in 1:B0) {            # loop begins for R0(X^*, F^hat)
  xyindex <- sample(seq(1:n), size=n, replace=T);
  tempx <- x[xyindex];
  tempy <- y[xyindex];
  temp <- lm(tempy~tempx)$coef;
  thetai <- temp[2]/temp[1];# theta^*_i
  R0star[i] <- thetai - theta;
  temp1 <- rep(0, B1);      # R0(X^**, F^hat)
  for(j in 1:B1) {          # loop begins for R1(X^*, F^hat)
    xyindex1 <- sample(seq(1:n), size=n, replace=T);
    tempx1 <- tempx[xyindex1];
    tempy1 <- tempy[xyindex1];
    temp <- lm(tempy1~tempx1)$coef;
    temp1[j] <- temp[2]/temp[1] - thetai;
  }
  R1star[i] <- sum(temp1 <= R0star[i])/B1;  # R1^hat(X^*, F^hat), (9.21) on page 267
}
proc.time() - tempt         # time cost
#   user  system elapsed 
# 517.05    0.34  543.31  
summary(R1star)
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
# 0.0000  0.3042  0.4917  0.5199  0.7567  1.0000
hist(R1star, xlab="R1(X^*, F^hat)", nclass=10)
xi1 <- quantile(R1star, c(0.025, 0.975))     # a 95% C.I. for R1^hat(X, F)
#     2.5%     97.5%
#0.0282500 0.9950833
temp <- quantile(R0star, xi1)                # a 95% C.I. for R0(X, F)
c(theta-temp[2], theta-temp[1])              # a 95% C.I. for theta=beta1/beta0
#         x          x 
#-0.1978562 -0.1664527