inst/thesis/fig/lgammafn_b.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

## ipeglim/inst/thesis/fig/lgammafn_b.R
## Chel Hee Lee <gnustats@gmail.com>
## Modified on 2013.NOV.07
##
## Note:
## my notation >> alpha=shape, beta=rate 
## R notation >> alpha=shape, beta=scale = 1/rate

rm(list=ls())

## Note:
## my notation >> alpha=shape, beta=rate >> E(Y) = alpha/beta
## R notation >> alpha=shape, beta=scale = 1/rate >> E(Y) = alpha*beta

## Left Panel: Log-Gamma Densities
alpha <- c(0.2, 0.5, 1, 2, 5)
mycolor <- rainbow(length(alpha))
beta <- 1/alpha

x <- seq(-5, 5, length.out=100)
y <- dlgamma(x, shape=1, rate=beta[1])

par(mfrow=c(1,2), mar=c(2,3,2,1)+0.1)
plot(x, y, type='l', xlim=c(-5, 5), ylim=c(0, 1), lwd=2,
     xlab=expression(theta), col=mycolor[1],
     ylab=expression(f(theta ~"|"~ alpha == 1 ~" ," ~ beta) ), 
     main="Log-Gamma")  
for(i in 2:5){
	curve(dlgamma(x, shape=1, rate=beta[i]), from=-5, to=5, add=TRUE, col=mycolor[i], lwd=2)
}
legend("topleft", legend=beta, lty=1, lwd=2, col=mycolor)

## Right panel: Gamma Densities
plot(x, dgamma(x, shape=1, rate=beta[1]), type="l", lwd=2,
     xlim=c(-5, 5), ylim=c(0, 1.2), col=mycolor[1],
     xlab=expression(theta), 
     ylab=expression(f(theta ~ "|" ~ alpha == 1 ~ ", " ~ beta)), 
     main="Gamma")
for(i in 2:5){
	curve(dgamma(x, shape=1, rate=beta[i]), from=-5, to=5, add=TRUE, col=mycolor[i], lwd=2)
}
legend("topleft", legend=beta, lty=1, lwd=2, col=mycolor)