Demo CLT (see how means get more concentrated, and their standard deviation shrinks to ~1/3 when n increases from 10 to 90)

print.r

# Number of samples to compute the mean. We try different n from 10, 20, ..., 90
n = seq(10, 90, 10)

# Number of means to compute the standard deviation of the mean
m = 30

set.seed(159)
# Generate a single random value from a "strange" distribution
gen.value <- function()
  sample(1:10, 1) +                               # Offset
  rbeta(1, 1, 10) * sample(1:10, 1) +             # Right-skewed distribution
  sample(c(rep(0, 50), 1:10)) * rnorm(1, 200, 10) # Outliers

all_ns = c()
all_means = c()
standard_deviations = c()

for (test_n in n) {
  # Mean generator
  gen.mean = function() mean(sapply(1:test_n, function(x) gen.value()))

  means = sapply(1:m, function(j) gen.mean())
  standard_deviations = append(standard_deviations, sd(means))
  all_means = append(all_means, means)
  all_ns = append(all_ns, rep(test_n, m))
}

pdf('print.pdf', width = 10, height = 5)
par(mfrow=c(1,2))
plot(all_means ~ all_ns)
plot(standard_deviations ~ n)
dev.off()