# Problem 4 of HW, STAT 789, Summer 2006

set.seed(321)
z <- rnorm(75, 0, 1)

# I'll generate 3200 bootstrap replicates of the sample median.  I'll use the first 800,
# the first 1600, and all 3200 of these to obtain 3 different bootstrap estimates of the
# standard error of the sample median, in order to check (approx) convergence.
samp.med.rep <- apply( matrix( sample( z, 75*3200, replace=T ), nrow=3200 ), 1, median )
nonpar.se.est <- numeric()
B <- c(800, 1600, 3200)
for ( b in 1:3 ) nonpar.se.est[b] <- sd( samp.med.rep[1:B[b]] )
part.a <- cbind( B, nonpar.se.est )
part.a
### I'll give 0.123 as the answer for part (a).

# I'll get the estimates of the parameters to use with nonparametric bootstrapping.
mean.est <- mean( z )
s.d.est <- sqrt( var( z )*74/75 )
mean.est; s.d.est
# Now I'll get the replicates.
samp.med.rep <- apply( matrix( rnorm( 75*3200, mean.est, s.d.est ), nrow=3200 ), 1, median )
par.se.est <- numeric()
for ( b in 1:3 ) par.se.est[b] <- sd( samp.med.rep[1:B[b]] )
part.b <- cbind( B, par.se.est )
part.b
### I'll give 0.141 as the answer for part (b).

# Now I'll get the replicates for the Monte Carlo approximation.
samp.med.rep <- apply( matrix( rnorm( 75*3200, 0, 1 ), nrow=3200 ), 1, median )
MC.se.est <- numeric()
for ( b in 1:3 ) MC.se.est[b] <- sd( samp.med.rep[1:B[b]] )
part.c <- cbind( B, MC.se.est )
part.c
### I'll give 0.142 as the answer for part (c).

asym.se <- sqrt( pi/150 )
asym.se
### I'll give 0.14472 as the answer for part (d).

# Now I'll paste in all of the code above that is after the generation of the
# original sample of size 75 to see how much the results depend of the random
# number seed used.