```# Goals: Scare the hell out of children with the Cauchy distribution.

# A function which simulates N draws from one of two distributions,
# and returns the mean obtained thusly.
one.simulation <- function(N=100, distribution="normal") {
if (distribution == "normal") {
x <- rnorm(N)
} else {
x <- rcauchy(N)
}
mean(x)
}

k1 <- density(replicate(1000, one.simulation(20)))
k2 <- density(replicate(1000, one.simulation(20, distribution="cauchy")))

xrange <- range(k1\$x, k2\$x)
plot(k1\$x, k1\$y, xlim=xrange, type="l", xlab="Estimated value", ylab="")
grid()
lines(k2\$x, k2\$y, col="red")
abline(v=.5)
legend(x="topleft", bty="n",
lty=c(1,1),
col=c("black", "red"),
legend=c("Mean of Normal", "Mean of Cauchy"))
# The distribution of the mean of normals collapses into a point;
# that of the cauchy does not.

# Here's more scary stuff --
for (i in 1:10) {
cat("Sigma of distribution of 1000 draws from mean of normal - ",
sd(replicate(1000, one.simulation(20))), "\n")
}
for (i in 1:10) {
cat("Sigma of distribution of 1000 draws from mean of cauchy - ",
sd(replicate(1000, one.simulation(20, distribution="cauchy"))), "\n")
}

# Exercise for the reader: Compare the distribution of the median of
# the Normal against the distribution of the median of the Cauchy.
```