# Bootstrap estimate of bias and standard deviation
# for tha_hat = max(x) where x[i] ~ U(0, tha)

set.seed(999)
tha <- 2
n <- 10
x <- runif(n, 0, tha)
cat('sample of size', n, 'from a U(0,tha), tha =', tha, '\n')

tha_hat <- max(x)
cat('tha estimate', tha_hat, '\n')

B <- 10000
thaB <- rep(0, B)
for (i in 1:B) {
  xB <- sample(x, n, replace=T)
  thaB[i] <- max(xB)
}
tha_bias <- mean(thaB) - tha_hat
cat('bootstrap estimate of bias', tha_bias, '\n')
tha_std <- sd(thaB)
cat('bootstrap estimate of std dev', tha_std, '\n')

N <- 10000
x <- matrix(runif(n*N, 0, tha), nrow=N, ncol=n)
tha_hat <- apply(x, 1, max)
tha_bias <- tha_hat - tha
cat('simulation estimate of bias', mean(tha_bias), '\n')
cat('simulation estimate of std dev', sd(tha_hat), '\n')