# These commands plot a centered gamma density
# and overlaid normal density with the same variance
par(mar=c(6,4,6,4),mfrow=c(2,2)) #Control margins, four plots/page
mu <- 2     # Set mean and dispersion of gamma 
phi <- 1/4  #
sh <- 1/phi # Shape and scale of gamma
sc <- mu/sh #
s <- sqrt(phi*mu^2)
y <- 6*s*(0:100)/100 # Set range of Y, zero to six sds.
x <- y-mu   # Transform to X=y-mu
plot(x,dgamma(x+mu,shape=sh,scale=sc),type="l",xlim=c(-4,4),ylim=c(0,0.5),ylab="Density")
text(-4,0.45,adj=0,expression(mu==2.00))
text(-4,0.40,adj=0,expression(phi==0.25))
x <- -4+(0:80)/10 # Set range of x from - 4 to 4
lines(x,dnorm(x,0,s),lty=2)
#
# To print the plot create a postscript file called "plot.eps"
postscript("plot.eps",horizontal=FALSE,onefile=FALSE)
par(mar=c(6,4,6,4),mfrow=c(2,2)) #Control margins, four plots/page
x <- y-mu   
plot(x,dgamma(x+mu,shape=sh,scale=sc),type="l",xlim=c(-4,4),ylim=c(0,0.5),ylab="Density")
text(-4,0.45,adj=0,expression(mu==2.00))
text(-4,0.40,adj=0,expression(phi==0.25))
x <- -4+(0:80)/10 
lines(x,dnorm(x,0,s),lty=2)
dev.off()
#You can view it (and then print it) by typing the command
#gv plot.eps &
#in an xterm window (after moving to the right directory).