## Simulate R random samples of size n from a standard normal population
n <- 41
R <- 10000
x <- matrix(rnorm(R * n), ncol = n)

## Compute sample means and sample medians
mn <- apply(x,1,mean)
md <- apply(x, 1, median)

## plot estimates of the spling distribution densities
plot(density(mn))
lines(density(md), col = "red")

## add the exaxt sampling distribution densities
z <- seq(-1,1,len=101)
lines(z, dnorm(z,sd=sqrt(1/n)), lty = 2)
stopifnot(n %% 2 == 1) # make sure n is odd
k <- (n - 1) / 2 + 1 
lines(z, dbeta(pnorm(z), k, n - k + 1) * dnorm(z), lty=2,col="red")

## estimate standard deviations of the sampling distributions
sd(mn)
sd(md)

## repeat in a loop for a range of n values
nn <- c(11, 21, 41, 61, 81, 101)
val <- matrix(0, nrow = length(nn), ncol = 2)
for (i in seq(along = nn)) {
    n <- nn[i]
    x <- matrix(rnorm(R * n), ncol = n)
    val[i, 1] <- sd(apply(x,1,mean))
    val[i, 2] <- sd(apply(x, 1, median))
}
colnames(val) <- c("mean", "median")

## some plots
matplot(nn, val, type = "l")
plot(nn, val[,2]/val[,1], type = "l")
plot(nn, val[,2]/val[,1], type = "l", ylim = c(0,2))
plot(nn, val[,2]/val[,1], type = "l", ylim = c(1,1.5))
plot(nn, val[,2], type="l")
lines(nn, 1.25/sqrt(nn), col = "red")