# First we compute the simulated results; reps=10,000:

results.med.est <- replicate(reps, med.est(sample(1:N,n)))
results.mean.est <- replicate(reps, mean.est(sample(1:N,n)))
results.max.est <- replicate(reps, max.est(sample(1:N,n)))
results.avgap2.est <- replicate(reps, avgap2.est(sample(1:N,n)))
results.iqr.est <- replicate(reps, iqr.est(sample(1:N,n)))
results.range.est <- replicate(reps, range.est(sample(1:N,n)))
results.norm.est <- replicate(reps, norm.est(sample(1:N,n)))

# Then we compute means of each data set:

m1 <- mean(results.med.est)
m2 <- mean(results.mean.est)
m3 <- mean(results.max.est)
m4 <- mean(results.avgap2.est)
m5 <- mean(results.iqr.est)
m6 <- mean(results.range.est)
m7 <- mean(results.norm.est)

# Next we compute the MSE of each data set:

mse1 <- mean((results.med.est - N)^2)
mse2 <- mean((results.mean.est -N)^2)
mse3 <- mean((results.max.est - N)^2)
mse4 <- mean((results.avgap2.est -N)^2)
mse5 <- mean((results.iqr.est - N)^2)
mse6 <- mean((results.range.est -N)^2)
mse7 <- mean((results.norm.est - N)^2)

# Now we form vectors of the means and MSEs, rounded:

means <- round(c(m1,m2,m3,m4,m5,m6,m7),2)
mse <- round(c(mse1,mse2,mse3,mse4,mse5,mse6,mse7),2)

# Finally, we ask R for a table:

cbind(means,mse)

estimator  means    mse

median     50.12   323.93
mean       50.18   152.78
max        49.80    65.66
avgp       49.92    69.76
IQR        50.82   247.67
Range      49.78   109.94
Norm       67.58   473.79



# Note the first column (estimator) I put in by hand.