Wednesday, January 18, 2012

Parameterizing a gamma distribution by mode and sd

When trying to fashion a gamma-shaped prior, I've found it more intuitive to start with the mode and standard deviation, instead of the mean and standard deviation as used in the book. The reason is that the gamma distribution is typically very skewed, and therefore the location of the mean is not very intuitive. Here I present how to compute the shape and rate parameters of the gamma distribution from a desired mode and standard deviation.

You can extract the mathematical formula for getting to shape and rate parameters from the R code I created:

# Specify desired mode and sd of gamma distribution:
mode = 10
sd = 20

# Here are the corresponding rate and shape parameter values:
ra = ( mode + sqrt( mode^2 + 4*sd^2 ) ) / ( 2 * sd^2 )
sh = 1 + mode * ra

show(sh)
show(ra)

# Graph it:
x = seq(0,mode+5*sd,len=1001)
plot( x , dgamma( x , shape=sh , rate=ra ) , type="l" ,
      main=paste("dgamma, mode=",mode,", sd=",sd,sep="") ,
      ylab=paste("dgamma( shape=",signif(sh,3)," , rate=",signif(ra,3)," )",
                 sep="") )
abline( v=mode , lty="dotted" )


The figure it creates looks like this (below). Notice that the mode is indeed as desired.



Hope that helps. Let me know.

3 comments:

  1. Often the gamma distribution is used as the prior for the precision of a normal distribution. One parametrization that I would like to do is to is specify the shape and rate using a mean SD and SD of the SD as if specifying an inverted gamma distribution as a prior for the SD of the normal distribution. It would be very convenient if this was possible because I feel it is easier to specify the variance prior on the same "scale" as the data, and not on the precision scale.

    I've tried to make this work, but I'm not sure I succeeded. My attempt is summarized in this question on Cross Validated. Any feedback and suggestions are greatly appreciated!

    /Rasmus

    ReplyDelete

  2. My impression (perhaps mistaken?) is that your goal is to put a prior on sigma instead of on tau (which equals 1/sigma^2), because it is more intuitive to deal with sigam. As Erik replied earlier on Cross Validated, this is straight forward in JAGS/BUGS:

    tau <- pow(sigma,-2)
    sigma ~ thePriorOfYourChoice

    There are various examples of this in Doing Bayesian Data Analysis, such as Figure 18.1, p. 494. For an example of putting a gamma prior on sigma, see this blog post: http://doingbayesiandataanalysis.blogspot.com/2012/04/improved-programs-for-hierarchical.html

    ReplyDelete
  3. I think I solved it now with the help of a mathematician friend. The solution was much more difficult than I could ever imagine. My solution can be found here.

    ReplyDelete