inst/thesis/fig/lgammafn_a.R

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

rm(list=ls())

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

x <- seq(-5, 5, length.out=100)
y <- dlgamma(x, shape=alpha[1], rate=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), col=mycolor[1],
     xlab=expression(theta), 
     ylab=expression(f(theta ~"|"~ alpha ~" ," ~ beta == 1) ),
     main="Log-Gamma")
for(i in 2:5){
	curve(dlgamma(x, shape=alpha[i], rate=1), from=-5, to=5, add=TRUE, col=mycolor[i], lwd=2)
}
legend("topleft", legend=alpha, lty=1, lwd=2, col=mycolor)

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

Try the ipeglim package in your browser

Any scripts or data that you put into this service are public.

ipeglim documentation built on May 2, 2019, 4:31 p.m.