## gamma-4.R 
## PJ Feb 16, 2002 / 2014-08-07
## This is a template showing how to see the central
## limit theorem as it applies to a gamma distribution

## You should alread have run template_gamma1 and 2.
## So you know the shape parameter is important.


## Run this a few times, changing the value of sh at
## the top. What do you see?

## Reminder: the mean of a gamma distribution is
## shape*scale.

##------------------------------------------##
## ---------RUN THIS PART FIRST ------------##
sh <- 0.2
sc <- 1
sampleSize <- 500
nOfRuns <- 2000

op <- par(no.readonly = TRUE)

SAVEME <- FALSE


getGammaMeans <- function (x){
  e <- rgamma(sampleSize, shape = sh, scale = sc)
  mean(e)
}

samplingDist <- sapply(1:nOfRuns, getGammaMeans)


if (SAVEME) pdf(file = "distributions-gamma-04-1.pdf", width = 8.5,
                height = 10, onefile = FALSE, family = "Times", paper = "special")


mytitle <- paste("The Sampling Distribution of a Gamma(sh=",sh,",sc=",sc,")")
hist(samplingDist, breaks=10, main = mytitle)

if (SAVEME) dev.off()

##------------------------------------------------##
# ---------- SECOND PART STARTS HERE -------------- #

# If you did not run that first part yet, go back.
# Then run this really fun part!

if (SAVEME) pdf(file = "distributions-gamma-04-2.pdf", width = 8.5,
                height = 10, onefile = FALSE, family = "Times", paper = "special")

## I've done it. It works.
par(mfrow = c(2,1))

oneGamma <- rgamma(sampleSize, shape = sh, scale = sc)

oneHist <- hist(oneGamma,
                main = paste("here is one sample from your Gamma(sh=",sh,",sc=",sc,")"))

text(.6*max(oneGamma), .8*max(oneHist$counts),
     paste("E(x)=", sc*sh, "\nobserved mean=",
           round(mean(oneGamma),3)), pos = 2)

mytitle <- paste("The Sampling Distribution of the Mean of Gamma(sh=",sh,",sc=",sc,")")
samplinghist <- hist(samplingDist,breaks=10,main=mytitle)

text(max(samplingDist), 0.8*max(samplinghist$counts),
     paste("The mean is", round(mean(samplingDist),3), "\nstd dev is ",
           round(sd(samplingDist),3)), pos=2)

if (SAVEME) dev.off()