# Problem 10 of HW,  STAT 789,  Summer 2006

set.seed(321)

usual.a <- numeric()
bs.a <- numeric()
jk.a <- numeric()
bs.b <- numeric()
jk.b <- numeric()
mean.jk.rep <- numeric()

for (j in 1:10000) {

  x <- rnorm(20, 1, 1)

  usual.a[j] <- sd(x)/sqrt(20)

  for (i in 1:20) mean.jk.rep[i] <- mean(x[-i])
  sq.mean.jk.rep <- mean.jk.rep^2
  jk.a[j] <- sd(mean.jk.rep)*19/sqrt(20)
  jk.b[j] <- sd(sq.mean.jk.rep)*19/sqrt(20)

  mean.bs.rep <- apply( matrix( sample( x, 20*800, replace=T ), nrow=800 ), 1, mean )
  sq.mean.bs.rep <- mean.bs.rep^2
  bs.a[j] <- sd(mean.bs.rep)
  bs.b[j] <- sd(sq.mean.bs.rep)

 }

# actual s.e. (a)
1/sqrt(20)

# mean usual estimates:
mean(usual.a)

# mean bootstrap estimates:
mean(bs.a)

# mean jackknife estimates:
mean(jk.a)

# s.d. usual estimates:
sd(usual.a)

# s.d. bootstrap estimates:
sd(bs.a)

# s.d. jackknife estimates:
sd(jk.a)

# est. MSE for bootstrap estimator:
bs.mse.a <- mean( ( bs.a - 1/sqrt(20) )^2 )
bs.mse.a

# est. MSE for jackknife estimator:
jk.mse.a <- mean( ( jk.a - 1/sqrt(20) )^2 )
jk.mse.a

## est. MSE (bootstrap) / est. MSE (jackknife)
bs.mse.a/jk.mse.a


# actual s.e. (b)
sqrt(0.205)

# mean bootstrap estimates:
mean(bs.b)

# mean jackknife estimates:
mean(jk.b)

# s.d. bootstrap estimates:
sd(bs.b)

# s.d. jackknife estimates:
sd(jk.b)

# est. MSE for bootstrap estimator:
bs.mse.b <- mean( ( bs.b - sqrt(0.205) )^2 )
bs.mse.b

# est. MSE for jackknife estimator:
jk.mse.b <- mean( ( jk.b - sqrt(0.205) )^2 )
jk.mse.b

## est. MSE (bootstrap) / est. MSE (jackknife)
bs.mse.b/jk.mse.b

###########################################################
###  Summary and Conclusions

## part (a)

# The estimand is 1/sqrt(20) = 0.2236

# The mean of the bootstrap estimates is 0.2149 and 
# the mean of the jackknife estimates is 0.2204.  So
# the bootstrap estimator seems to have a larger bias.
# Whereas the jackknife estimate for a given sample is 
# s/sqrt(n), which is the usual estimate of the standard 
# error of the sample mean, the bootstrap estimate is 
# closer to s*/sqrt(n), where s* denotes the plug-in 
# estimate of the distribution standard deviation.
# The plug-in estimate, with n in the denominator 
# instead of n-1, isn't as good.

# The standard deviation of the bootstrap estimates is 
# 0.0352 and the standard deviation of the jackknife 
# estimates is 0.0357.  So there is little difference
# in the standard deviations.  

# While for each estimator, its standard error is 
# appreciably larger than its bias, when comparing the
# two estimators, the difference in bias is much
# larger than the difference in the standard errors.
# The estimated MSE for the bootstrap estimator is
# more than 2% larger than the estimated MSE for the
# jackknife estimator.

# Since the sample mean is a linear statistic, one
# expects the jackknife estimator to be competitive with
# the bootstrap estimator, and in this case it's better
# because of the bootstrap estimator's greater bias.

## part (b)

# The estimand is sqrt(0.205) = 0.4528.

# The mean of the bootstrap estimates is 0.4344 and 
# the mean of the jackknife estimates is 0.4400.  So
# the bootstrap estimator seems to have a larger bias.

# The standard deviation of the bootstrap estimates is 
# 0.1224 and the standard deviation of the jackknife 
# estimates is 0.1243.  So the jackknife estimator
# seems to have a larger standard error.

# Since the square of the sample mean is a quadratic 
# statistic, one expects that the bootstrap estimate
# of standard error will be superior to the jackknife
# estimate of standard error.  When comparing estimated
# MSEs, we find that the bootstrap estimator is about 2% 
# better.  (The jackknife estimator of standard error did 
# better with regard to bias, and the bootstrap estimator 
# of standard error did better with regard to standard error.
# But overall, when comparing the estimated MSEs of the two
# estimators, the bootstrap estimator is the winner.)