###################################################
### chunk number 1: startup
###################################################
    rm(list=ls())
    set.seed(20090804)
    options(width=55,digits=4,continue=" ")
    par(ask=FALSE, mfrow=c(1,1))


###################################################
### chunk number 2: aircond
###################################################
  library(boot)
  t(aircondit7)     #7th aircraft


###################################################
### chunk number 3: cities
###################################################
  t(head(bigcity, n=15))


###################################################
### chunk number 4: bcityplot
###################################################
    plot(bigcity,
      xlab="1920 pop.", ylab="1930 pop.", cex=.6)


###################################################
### chunk number 5: 
###################################################
library(boot)
y <- aircondit7$hours
n <- length(y)
ybar <- mean(y)


###################################################
### chunk number 6: qqplot
###################################################
  q <- qexp(ppoints(24), rate=1/ybar)
  qqplot(y, q)
  abline(0, 1)


###################################################
### chunk number 7: 
###################################################
 R <- 2000
 Tr <- replicate(R, expr={
            x <- rexp(n, rate=1/ybar)
            mean(x)})
 bias <- mean(Tr) - ybar
 bias


###################################################
### chunk number 8: 
###################################################
 se <- sd(Tr)
 data.frame(bias=bias, se=se, var=se^2)


###################################################
### chunk number 9: 
###################################################
vy <- ybar^2 / n
sqrt(vy)


###################################################
### chunk number 10: ftimeplot
###################################################
    y <- aircondit7$hours
    plot(ecdf(y))
    curve(pexp(x, rate=1/mean(y)), add=TRUE)


###################################################
### chunk number 11: varMLE
###################################################
    s.hat <- function(x) {
        m <- mean(x)
        sqrt(mean((x - m)^2))
    }


###################################################
### chunk number 12: 
###################################################
 R <- 2000
 ybar <- mean(y)
 Tr <- replicate(R, expr={
            x <- sample(y, size=n, replace=TRUE)
            mean(x)})
 bias.boot <- mean(Tr) - ybar
 se.boot <- s.hat(Tr)
 est <- data.frame(bias=bias.boot, se=se.boot)
 rbind(est, c(bias, se), c(0, sqrt(vy)))


###################################################
### chunk number 13: ratio
###################################################
    ratio <- function(X) mean(X[,2])/mean(X[,1])


###################################################
### chunk number 14: 
###################################################
    data(city)
    y <- as.matrix(city)
    t(y)
    i <- sample(1:10, replace=TRUE)
    x <- y[i, ]
    t(x)


###################################################
### chunk number 16: cities
###################################################
    data(bigcity)
    y <- as.matrix(bigcity)
    T0 <- ratio(y)
    R <- 999
    Tr <- replicate(R, expr = {
        i <- sample(1:NROW(y), replace=TRUE)
        x <- y[i, ]
        ratio(x)})
    c(T0, mean(Tr), s.hat(Tr))


###################################################
### chunk number 17: ratiohist
###################################################
    hist(Tr, prob=TRUE, breaks="scott",
        main="Replicates of ratio, R=999")
    lines(density(Tr))
    points(T0, 0, pch=17, cex=2)


###################################################
### chunk number 18: ratio.boot
###################################################
    ratio.boot <- function(dat, i) {
        x <- dat[i, ]
        mean(x[, 2]) / mean(x[, 1])
    }


###################################################
### chunk number 19: boot-ratio
###################################################
    # library(boot)   # if not loaded
    y <- as.matrix(bigcity)
    boot(y, statistic=ratio.boot, R=999)


###################################################
### chunk number 20: boot-object
###################################################
boot.obj <- boot(y, statistic=ratio.boot, R=999)
names(boot.obj)
boot.obj$t0                  #observed statistic T
head(as.vector(boot.obj$t))  #replicates of T


###################################################
### chunk number 21: quantiles
###################################################
    quantile(boot.obj$t, c(.05, .95), type=1)


###################################################
### chunk number 22: normal
###################################################
alpha <- c(.025, .975)
boot.obj$t0 + qnorm(alpha) * sd(boot.obj$t)


###################################################
### chunk number 23: percentile
###################################################
quantile(boot.obj$t, alpha, type=6)


###################################################
### chunk number 24: basic-bootstrap
###################################################
2*boot.obj$t0 -
     quantile(boot.obj$t, rev(alpha), type=1)


###################################################
### chunk number 25: bootci
###################################################
boot.ci(boot.obj, type=c("norm","perc","basic"))


###################################################
### chunk number 26: bootci
###################################################
ci <- boot.ci(boot.obj, type=c("norm","perc","basic"))
names(ci)
ci$percent
L <- ci$percent[4]
U <- ci$percent[5]
c(L, U)


###################################################
### chunk number 27: 
###################################################
y <- as.matrix(bigcity)
boot.out <- boot(y, ratio.boot, R=999)
ci <- boot.ci(boot.out, type=c("perc","bca"))
ci