n<-1e6
p<-5

keep<-function(x,rng){
x[(rng[1]<x)&(x<rng[2])]
}

pdf(file="cauchy.xbar.pdf")

par(mfrow=c(2,1))

# Normal
x<-rnorm(n)
xbar<-rowMeans(matrix(rnorm(n*p),n,p))

xlim<-c(-3,3)
x<-keep(x,xlim)
xbar<-keep(xbar,xlim)

breaks<-seq(-3,3,length=51)
hist.x<-hist(x,plot=FALSE,breaks=breaks)
hist.xbar<-hist(xbar,plot=FALSE,breaks=breaks)

ylim<-c(0,max(hist.x$density,hist.xbar$density))
plot(hist.x,xlim=xlim,ylim=ylim,freq=FALSE,
     main="Histogram of X and Xbar for Normal Population",xlab="")
plot(hist.xbar,xlim=xlim,ylim=ylim,freq=FALSE,border="red",add=TRUE)

# Cauchy
x<-rcauchy(n)
xbar<-rowMeans(matrix(rcauchy(n*p),n,p))

xlim<-c(-3,3)
x<-keep(x,xlim)
xbar<-keep(xbar,xlim)

breaks<-seq(-3,3,length=51)
hist.x<-hist(x,plot=FALSE,breaks=breaks)
hist.xbar<-hist(xbar,plot=FALSE,breaks=breaks)

ylim<-c(0,max(hist.x$density,hist.xbar$density))
plot(hist.x,xlim=xlim,ylim=ylim,freq=FALSE,
     main="Histogram of X and Xbar for Cauchy Population",xlab="")
plot(hist.xbar,xlim=xlim,ylim=ylim,freq=FALSE,border="red",add=TRUE)

dev.off()