Suppose we want to simulate for , . Suppose we want to do the same for the Cauchy distribution.

In other words, we want to draw several random variables from a normal distribution and then take the average. As increases we should get closer to the mean of the distribution we’re drawing from, 0 in this case.

The R code below will do this. It produces this graph:

Notice that while the Normal distribution converges quickly the Cauchy never does. This is because the Cauchy distribution has fat tails and so extreme observations are common.

################################################################
# R Simulation
# James McCammon
# 2/20/2017
################################################################
# This script goes through the simulation of plotting both normal
# and Cauchy means for random vectors of size 1 to 10,000. It
# also demonstrates function creation and plotting.
# Highlight desired section and click "Run."
# Set working directory as needed
setwd("~/R Projects/All of Statistics")
###
# Calculate means and plot using base R
###
# Set seed for reproducibility
set.seed(271)
#
# Version 1: Simple
#
n = seq(from=1, to=10000, by=1)
y=sapply(n, FUN=function(x) sum(rnorm(x))/x)
plot(n, y, type="l")
#
# Version 2: Define a function
#
sim = function(x, FUN) {
sapply(x, FUN=function(x) sum(FUN(x))/x)
}
# Use function to plot normal means
x = seq(from=1, to=10000, by=1)
y1 = sim(x, rnorm)
plot(x, y1, type="l")
# Use function to plot Cauchy means
y2 = sim(x, rcauchy)
plot(x, y2, type="l")
#
# Version 3: More complex function
#
# This function has:
# (1) error checking
# (2) extra argument options
# (3) the ability to input any distribution R supports
sim = function(x, FUN, ...) {
if(!is.character(FUN)) stop('Please enter distribution as string.')
dists = c('rnorm',
'rbeta',
'rbinom',
'rcauchy',
'rchisq',
'rexp',
'rf',
'rgamma',
'rgeom',
'rhyper',
'rlnorm',
'rmultinom',
'rnbinom',
'rpois',
'rt',
'runif',
'rweibull')
if(is.na(match(FUN, dists))) stop(paste('Please enter a valid distribution from one of:', paste(dists, collapse=', ')))
FUN = get(FUN)
sapply(x, FUN=function(x) sum(FUN(x, ...))/x)
}
# We have to define our function in string form.
# This will throw error 1.
test1 = sim(x, rnorm)
# We have to input a distribution R supports.
# This will throw error 2.
test2 = sim(x, 'my_cool_function')
# We can input additional arguments like the
# mean, standard deviations, or other shape parameters.
test3 = sim(x, 'rnorm', mean=10, sd=2)
####
# Using ggplot2 to make pretty graph
###
# Load libraries
library(ggplot2)
library(ggthemes)
library(gridExtra)
png(filename='Ch3-Pr9.png', width=1200, height=600)
df1 = cbind.data.frame(x, y1)
p1 = ggplot(df1, aes(x=x, y=y1)) +
geom_line(size = 1, color='#937EBF') +
theme_fivethirtyeight() +
ggtitle("Normal Means")
df2 = cbind.data.frame(x, y2)
p2 = ggplot(df2, aes(x=x, y=y2)) +
geom_line(size = 1, color='#EF4664') +
theme_fivethirtyeight() +
ggtitle("Cauchy Means")
# Save charts
grid.arrange(p1,p2,nrow=2,ncol=1)
dev.off()
# Print charts to screen
grid.arrange(p1,p2,nrow=2,ncol=1)