## Problem 1

reps <- 1e4

for (n in c(5,10,20))
{
  X <- array(rnorm(n*reps,1,1), c(n,reps))
  sampMeans <- apply(X, 2, mean)
  sampMeanSq <- mean((sampMeans)^2)
  print(sampMeanSq)
}


## Problem 2
reps <- 1e4

nresamp <- 1000

X <- rnorm(n)
p <- length(X)
ii <- ceiling(p*runif(p*nresamp)) 
B <- X[ii]
B <- array(B, c(p,nresamp))




# Problem 3
n <- 20
X <- rnorm(n,1,1)
xbar <- mean(X)
sigmahat <- sd(X)

## Bootstrap samples.
ii <- ceiling(n*runif(1000*n))
B <- array(X[ii], c(n,1000))
  
xbar_k <- apply(B, 2, mean)
sigmahat_k <- apply(B, 2, sd)
MD <- abs(xbar_k-xbar) / sigmahat_k
MDsorted <- sort(MD)
F <- MDsorted[950]

c1 <- xbar - F*sigmahat
c2 <- xbar + F*sigmahat
print(c(c1,c2))