###### Example of Nonparametric Bootstrap for City Population Data ######
#
# For use with city data from Davison and Hinkley
#

library(boot)
attach(city)
city


nonparboot.fn <- function(B, alpha){

index <- c(1:10)

n <- nrow(city)

M <- matrix(data = NA, nrow = nrow(city), ncol = ncol(city)) 
t <- mean(x)/mean(u)

tstar <- c(rep(NA, B))

for(j in 1:B){

	### Nonparametric Bootstrap ###	

	for(k in 1:n){

		i <- sample(index, size = 1, replace = TRUE)
		M[k,1] <- city[i,1]
		M[k,2] <- city[i,2] 
	
	}

	tstar[j] <- mean(M[,2])/mean(M[,1])

}

var.delta <- sum((x - t*u)^2)/(n*mean(u))^2

mean <- mean(tstar)
bias <- mean(tstar) - t
variance <- var(tstar)

### ordered tstar values ###
ts <- sort(tstar)

### 1-2*alpha equal-tail CI (2.10) in Davison & Hinkley ###

lower <- t - (ts[(B+1)*(1-alpha)] - t) 
upper <- t - (ts[(B+1)*alpha] - t)

print(cbind(tstar))

print(cbind(mean, bias, variance, var.delta))
print(cbind(lower, upper))


hist(tstar, freq = FALSE, main = "Bootstrap Distribution")

}

nonparboot.fn(999, 0.025)