### ASSIGNMENT 1, PROBLEM 3
### 
### In R you can do
### 
### source("gamma-normal-approx.R")
### 
### to see the results of this script.
### 
### For printing, use the commented out commands at the bottom to
### produce a pdf file.
### 

dgamma2 <-
  function(y, phi, mu=1, ...) dgamma(y, shape = 1/phi, scale = phi * mu, ...)
qgamma2 <-
  function(p, phi, mu=1, ...) qgamma(p, shape = 1/phi, scale = phi * mu, ...)

fungamplot <- function(fun, gfuninv,
                       zlim = 4, npt = 100, tail.omit = pnorm(-abs(zlim)),
                       log.grid = FALSE, maxfactor = 4) {
  opar <- par(mfrow=c(3,3))
  zlim <- abs(zlim)
  for (i in -1:7) {
    alpha <- 2^i
    phi <- 1/alpha
    ymin <- qgamma2(tail.omit, phi)
    ymax <- qgamma2(1 - tail.omit, phi)
    if (log.grid) y <- exp(seq(log(ymin), log(ymax), length=npt))
    else y <- seq(ymin, ymax, length=npt)
    x <- fun(y, phi)
    fx <- dgamma2(y, phi) * abs(gfuninv(y, phi))
    z <- seq(-zlim, zlim, length=npt)
    fz <- dnorm(z)
    fmax <- max(fx, fz)
    fmax <- min(fmax, maxfactor*dnorm(0))
    plot(c(-zlim,zlim), c(0, fmax), xlab="x", ylab="f(x)", type="n")
    lines(x, fx)
    lines(z, fz, lty=2)
    mtext(bquote(phi == 2^.(-i)))
  }
  par(opar)
  NULL
}

fun3b <- function(y, phi) (y - 1)/sqrt(phi)

gfuninv3b <- function(y, phi) sqrt(phi)

fungamplot(fun3b, gfuninv3b)

## 
## To produce a nice pdf version, try these commands
## 
##   pdf("asgn1-p3b-plot.pdf", onefile=FALSE, width=6, height=6, pointsize=10)
##   par(mar=c(4,4,2,1)+0.1, mgp=c(2,0.75,0), oma=c(0,0,0,0)+0.1)
##   fungamplot(fun3b, gfuninv3b, npt=100)
##   dev.off()
## 
## For part (d) you will want to use a call like
##
##    fungamplot(fun3d, gfuninv3d, log.grid = TRUE)
##
## but you have to write fun3d and gfuninv3d first.  Ditto for part
## (e).  You may want to try setting log.grid to false here to see why
## I added it to the function --- feel free to change it if you come
## up with something better.