# STATS 4CI3 Winter 2019

# Example of smoothed and shrunk-smoothed bootstrap

# Let us consider the median of the airconditioning data
# First we will examine the usual nonparametric bootstrap

med <- median(times)
R <- 10000
n <- length(times)
times.boot1 <- matrix(sample(times, n*R, replace=T), ncol=n)
med.star1 <- apply(times.boot1, 1, median)

Ystar <- matrix(rnorm(R*n), ncol=n)
h <- density(times)$bw
times.boot2 <- times.boot1+h*Ystar
med.star2 <- apply(times.boot2, 1, median)

b <- 1/sqrt(1+n*h^2/((n-1)*var(times)))
a <- (1-b)*mean(times)
times.boot3 <- a+b*times.boot1+h*b*Ystar
med.star3 <- apply(times.boot3, 1, median)

colMeans(cbind(med.star1, med.star2, med.star3))-med
library(matrixStats)
colSds(cbind(med.star1, med.star2, med.star3))

######################################################

# Bootstrap Confidence Intervals
lambda.hat <- 1/mean(times)

R <- 10000
n <- length(times)

xstar <- matrix(sample(times,n*R, replace=TRUE), ncol=n)
lambda.hat.star <- 1/rowMeans(xstar)

b.boot <- mean(lambda.hat.star-lambda.hat)
v.boot <- var(lambda.hat.star)

# Bootstrap Normal Confidence Interval
lambda.hat-b.boot+c(-1,1)*qnorm(0.975)*sqrt(v.boot)

# Basic Bootstrap Confidence Interval
2*lambda.hat-sort(lambda.hat.star)[c(0.975,0.025)*R]

# Work on the log scale
phi.hat <- log (lambda.hat)
phi.hat.star <- log(lambda.hat.star)

b.boot.log <- mean(phi.hat.star-phi.hat)
v.boot.log <- var(phi.hat.star)

# A 95% bootstrap normal CI for phi=log(lambda) is then
phi.hat-b.boot.log+c(-1,1)*qnorm(0.975)*sqrt(v.boot.log)

# Transforming back to the lambda scale we get
exp(phi.hat-b.boot.log+c(-1,1)*qnorm(0.975)*sqrt(v.boot.log))

# Similarly the basic bootstrap intervals
2*phi.hat-sort(phi.hat.star)[c(0.975,0.025)*R]
exp(2*phi.hat-sort(phi.hat.star)[c(0.975,0.025)*R])

# Percentile Interval
sort(lambda.hat.star)[c(0.025,0.975)*R]

###############################################################