# To reproduce the wikipedia page graphic
x <- seq(from=0, to=8, by=0.1)
plot(x, dggamma(x, theta=2, kappa=0.5, delta=0.5), lty=1, col="blue",
type="l", lwd=2, xlab="x", ylab="PDF")
lines(x, dggamma(x, theta=1, kappa=1, delta=0.5), lty=1, col="green", lwd=2)
lines(x, dggamma(x, theta=2, kappa=1, delta=2), lty=1, col="red", lwd=2)
lines(x, dggamma(x, theta=5, kappa=1, delta=5), lty=1, col="yellow", lwd=2)
lines(x, dggamma(x, theta=7, kappa=1, delta=7), lty=1, col="grey", lwd=2)
legend("topright", legend=c("a=2, d=0.5, p=0.5", "a=1, d=1, p=0.5",
"a=2, d=1, p=2", "a=5, d=1, p=5", "a=7, d=1, p=7"),
col=c("blue", "green", "red", "yellow", "grey"),
lty=1, lwd=2, bty="n")
par <- c(theta=2, kappa=0.5, delta=0.5)
# Mean, var and sd
mean.ggamma <- function(theta, kappa, delta)
return(theta*(gamma((kappa+1)/delta))/gamma(kappa/delta))
var.ggamma <- function(theta, kappa, delta)
return(theta^2* ( ( (gamma((kappa+2)/delta))/gamma(kappa/delta) ) -
( (gamma((kappa+1)/delta))/gamma(kappa/delta) )^2 ) )
sd.ggamma <- function(theta, kappa, delta)
return(sqrt(theta^2* ( ( (gamma((kappa+2)/delta))/gamma(kappa/delta) ) -
( (gamma((kappa+1)/delta))/gamma(kappa/delta) )^2 ) ))
Run the code above in your browser using DataLab