Created
October 7, 2015 15:19
Fit a gamma distribution based on prior information
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment