###########################
#### bootstrap package ####
###########################

# package bootstrap

install.packages("bootstrap")
library(bootstrap)

############################
# jackknife(x, theta, ...) #
############################

# Example 1
set.seed(2014)
n=200
x<-rnorm(n)
theta<-function(x){mean(x)}
results<-jackknife(x,theta)
(mean(results$jack.values)-mean(x))*(n-1);results$jack.bias   # bias
sd(results$jack.values)*(n-1)/sqrt(n);results$jack.se         # se of theta

# Example 2
# using my own code
data(patch)
n<-nrow(patch)
y<-patch$y
z<-patch$z
theta.hat<-mean(y)/mean(z)
theta.jack<-numeric(n)
for(i in 1:n){theta.jack[i]<-mean(y[-i])/mean(z[-i])}
bias<-(n-1)*(mean(theta.jack)-theta.hat)
sd(theta.jack)*(n-1)/sqrt(n)
bias

# using the function
xdata<-cbind(y,z)
theta<-function(x,xdata){mean(xdata[x,1])/mean(xdata[x,2])}
results<-jackknife(1:n,theta,xdata)
results$jack.se
results$jack.bias

###########################################
# bootstrap(x,nboot,theta,..., func=NULL) #
###########################################

n<-20
B<-20000
set.seed(2014)
x<-rnorm(n)
theta<-function(x){mean(x)}      # gives us sample mean for each bootstrap sample
results<-bootstrap(x,B,theta)
ssdd<-function(x){sd(x)}         # gives us sd for the bootstrap sample mean
results<-bootstrap(x,B,theta,func=ssdd)
results$func.thetastar;sd(results$thetastar)   # sd of sample mean
results$jack.boot.se;sd(results$jack.boot.val)*(n-1)/sqrt(n)
#  jackknife-after-bootstrap standard error estimate

# If I did the jackknife-after-bootstrap standard error using my own code,
# it turns out to be very close to the results using the function

theta=numeric(B)
indices<-matrix(0,nrow=B,ncol=n)
for(b in 1:B){
              ind<-sample(1:n,size=n,replace=TRUE)   
              theta[b]=mean(x[ind]) 
              indices[b,]<-ind
              }
se.jack<-numeric(n)
for(i in 1:n){
              keep<-(1:B)[apply(indices,MARGIN=1,FUN=function(k){!any(k==i)})]
              se.jack[i]<-sd(theta[keep]) 
              }
sd(theta)
sqrt((n-1)*mean((se.jack-mean(se.jack))^2))

# we can also use this to estimate the 95th percentile of the sample mean 
# and its jackknife-after-bootstrap standard error

perc95 <- function(x){quantile(x, .95)}
results <- bootstrap(x,100,theta, func=perc95)
results

###########################################################
# crossval(x, y, theta.fit, theta.predict, ..., ngroup=n) #
###########################################################

set.seed(2014)
x<-rnorm(85)
y<-2*x+.5*rnorm(85)
theta.fit<-function(x,y){lsfit(x,y)}
theta.predict<-function(fit,x){
                               cbind(1,x)%*%fit$coef
                               }
results<-crossval(x,y,theta.fit,theta.predict,ngroup=85)
results