### 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.