## Title: distributions-GammaVersusNormal-1.R
## Author: Paul Johnson <pauljohn at ku.edu>
## Date posted: 2013-02-08
## Description. Explore the Gamma PDF and its comparison
## to the Normal distribution. This example shows some
## nice uses of lines and legends in the plot at the end.


## Change to TRUE if you want to write the file
SAVEME <- FALSE

pdf.options(family = "Times", paper = "special", onefile = TRUE)

## Ask for a big sheet of paper to draw on:
if (SAVEME) pdf(file = "gammaVnormal-1.1.pdf", width = 8.5, height = 11)
## If you are running interactively, I suggest you make the
## window larger.

## Partition display into 5 rows, 2 columns
par(mfcol=(c(5,2)))

#first, just leave the gamma "shape" equal to 1.  Then the
# mean = scale and the variance= scale*scale
for (i in 1:10) {
    sc <- 0.5 + (i-1)
    title <- paste("Histogram, shape=1, scale=", sc)
    x <- rgamma(1250, shape = 1, scale = sc)
    hist(x, breaks = 60, main = title)
}

if (SAVEME) dev.off()

#what if we fix the scale=5 and shape varies?

if (SAVEME) pdf(file = "gammaVnormal-1.2.pdf", width = 8.5, height = 11)
par(mfcol=(c(5,2)))
for (i in 1:10) {
    sc <- 5
    sh = 0.5 + (i-1)
    title <- paste("Histogram, shape=,",sh,"scale=5")
    x <- rgamma(1250, shape = sh, scale = sc)
    hist(x, breaks = 60, main = title)
}

if (SAVEME) dev.off()



### Create a gamma and normal sample with same mean/SD
if (SAVEME) pdf(file = "gammaVnormal-1.3.pdf", width = 7, height = 7)
par(mfcol=c(1,1))
y <- rgamma(1000, scale = 10, shape = 1);
x <- rnorm(1000, mean = 10, sd = 10)
plot(x, y, main = "y a gamma,x a normal")

if (SAVEME) dev.off()

### here's a fun exercise. Compare the observed density
if (SAVEME) pdf(file = "gammaVnormal-1.4.pdf", width = 8.5, height = 11)
plot(density(y), ylim=c(0,.1), lty = 2, xlab = "x", main = "Gamma and Normal Density")
dennorm <- density(x)
lines(dennorm, lty = 3)
legend("topright", legend=c("E[x]=10,SD[x]=10",
                       "gamma sample density",
                       "normal sample density"), lty = c(NA,2,3))

if (SAVEME) dev.off()

### Now, start over, add 2 "true" density lines
if (SAVEME) pdf(file = "gammaVnormal-1.5.pdf", width = 7, height = 7)

plot(density(y), ylim = c(0,.1), lty = 1, xlab = "x",
     main = "Gamma and Normal Density")

dennorm <- density(x)
lines(dennorm, lty = 3)

xrange <- seq(range(y)[1], range(y)[2], length.out = 100)

truedengamma <- dgamma(xrange, scale = 10, shape = 1)
lines(xrange, truedengamma, lty = 4,col = "gray90", lwd = 5)

truedennorm <- dnorm(xrange, mean = 10, sd = 10)
lines(xrange, truedennorm, lty = 5, col = "lightpink", lwd = 5)

legend("topright", legend = c("E[x]=10,SD[x]=10",
                       "gamma sample density", "normal sample density", "theoretical gamma",
                       "theoretical normal"), lty = c(NA,1,3,4,5), col =
       c(1,1,1,"gray90","lightpink"), lwd = c(NA,1,1,5,5))

if (SAVEME) dev.off()