# 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.