Fit a gamma distribution based on prior information

gamma_parameters.R

# Objective
#
# Fit a gamma distribution knowing that:
# - 20% fall below 15 days
# - 80% fall below 60 days
# Inspired by http://www.johndcook.com/blog/2010/01/31/parameters-from-percentiles/

x <- c(0.2, 0.8)
y <- c(15, 60)

# ?dgamma
# ?qgamma

# plot the gamma curve for a given shape and rate argument ----------------

# draws vertical lines at p
plotGamma <- function(shape=2, rate=0.5, to=0.99, p=c(0.1, 0.9), cex=1, ...){
  to <- qgamma(p=to, shape=shape, rate=rate)
  curve(dgamma(x, shape, rate), from=0, to=to, n=500, type="l", 
        main=sprintf("gamma(x, shape=%1.2f, rate=%1.2f)", shape, rate),
        bty="n", xaxs="i", yaxs="i", col="blue", xlab="", ylab="", 
        las=1, lwd=2, cex=cex, cex.axis=cex, cex.main=cex, ...)
  gx <- qgamma(p=p,  shape=shape, rate=rate)
  gy <- dgamma(x=gx, shape=shape, rate=rate)
  for(i in seq_along(p)) { lines(x=rep(gx[i], 2), y=c(0, gy[i]), col="blue") }
  for(i in seq_along(p)) { text(x=gx[i], 0, p[i], adj=c(1.1, -0.2), cex=cex) }
}

# plot several gamma curves -----------------------------------------------

# note the shape argument determines the overall shape
# the scale parameter only affects the scale values
oldpar <- par(mfrow = c(2, 3), mai=rep(0.5, 4))
plotGamma(1, 0.1)
plotGamma(2, 0.1)
plotGamma(6, 0.1)
plotGamma(1, 1)
plotGamma(2, 1)
plotGamma(6, 1)
par(oldpar)


# formulate a gamma error function for non-linear least squares -----------

errorGamma <- function(p=c(10, 3), quantiles, exp){
  gx <- qgamma(p=quantiles, shape=p[1], rate=p[2])
  sum((gx-exp)^2)
}

p <- c(0.1, 0.8)
solution <- nlm(f=errorGamma, p=c(2, 1), quantiles=p, exp=c(30, 90))

shape <- solution$estimate[1]
rate  <- solution$estimate[2]
plotGamma(shape, rate, p=p)