7
$\begingroup$

The gamma distribution is a commonly used prior distribution for the precision ($1/sd^2$) of a normal distribution in Bayesian hierarchical modeling. I want to use an informed prior for the variance and I ask an expert on the data I'm trying to model. An easy way of describing earlier knowledge about the variance of the data is using the standard deviation as it is on the same scale as the data. My expert tells me: "Well in general the data tend to have a SD of 100, and I believe the SD of the SD is at most 50". I believe it would be much harder for my expert directly give this statement using precision and the SD of the precision.

Now that I have the expert's statement I want to incorporate it into my JAGS/BUGS model but in JAGS/BUGS I can't specify the gamma prior directly using this data as gamma is parametrized using shape and rate and the normal distribution is parameterized using precision.

What I would want to do is to take the statements of the mean of the sampling distribution's SD ($\mu_\sigma$) and the SD of the sampling distribution's SD ($\sigma_\sigma$) and use these parameters to set up the corresponding gamma prior for the precision of the sampling distribution. How could I do that? That is, given the expert's guestimate of the mean SD and the SD of the SD of the data, how do I specify the corresponding hyperparameters shape and rate for the gamma prior on the sampling distributions precision?

The following diagram shows the model and illustrates my question:

enter image description here

The code below would be the corresponding BUGS/JAGS model, where the "functions" f1 and f2 are the functions I would like to know how to implement.

model {
  for(i in 1:length(y)) {
    y ~ dnorm(mu, tau)
  }

  tau ~ dgamma(shape, rate)
  shape <- f1(mu_sigma, sigma_sigma)
  rate <- f2(mu_sigma, sigma_sigma)
  mu_sigma <- 100
  sigma_sigma <- 50

  mu ~ dnorm(M, P)
  M <- 0
  P <- 0.0001
}

The reason why I would want to specify the shape and the rate of the gamma using $\mu_\sigma$ and $\sigma_\sigma$ is because I would find it more intuitive to think in the scale of the SD of the sampling distribution rather than in the scale of the precision.

$\endgroup$
  • $\begingroup$ I believe you can stipulate a prior for the SD directly. Just use the reciprocal of that SD's square when parameterizing the normal distribution. $\endgroup$ – whuber Oct 25 '12 at 14:07
  • $\begingroup$ Yes, but if I want to use a gamma distribution how do I do then? An example or a piece of JAGS/BUGS code would be much appreciated! $\endgroup$ – Rasmus Bååth Oct 25 '12 at 14:52
  • 1
    $\begingroup$ If you want to stipulate a prior for the SD do it like this: sd ~ someprior and prec <- pow(1/sd,2). Then use prec as the precision for your normal. $\endgroup$ – Erik Nov 6 '12 at 15:53
  • $\begingroup$ You would want to be careful and not make your prior too informative. In most cases, you would want very small values of both $\alpha$ and $\beta$ to ensure a reasonable center of the distribution (of which you may have some preliminary idea) yet allowing for the data to modify your belief in case it was wrong, and there is way more or way less variability than you guessed initially. $\endgroup$ – StasK Nov 6 '12 at 17:24
6
+100
$\begingroup$

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, 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

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ We appreciate your enthusiasm for the site, @JohnKruschke, & help w/ questions here. Please maintain 1 account, though. This will also help you, as you will be notified when someone leaves a comment, eg, to one of your posts. $\endgroup$ – gung - Reinstate Monica Nov 7 '12 at 7:25
  • 1
    $\begingroup$ Well, what I want to do is to use a gamma prior on tau. But instead of specifying the shape and rate of the gamma prior by specifying the mean and SD of tau I would like to specify the mean and SD of sigma and the calculate, in some way, the corresponding shape and rate for the gamma prior on precision. Another way would be (I guess) to find a "thePriorOfYourChoice" that gives tau a gamma distribution but that can be parametrized with the mean and SD of sigma. If this is difficult I would be happy for any suggestion for an informative variance prior parametrizable with mean and SD of sigma. $\endgroup$ – Rasmus Bååth Nov 7 '12 at 8:34
  • 1
    $\begingroup$ John, I've merged your 4 accounts (thanks to @gung who noticed that). Please, register your account once and for all so that you get all benefits from Stack Exchange system-wide facilities. $\endgroup$ – chl Nov 7 '12 at 9:49
6
$\begingroup$

I eventually worked out the answer myself (with the help of a mathematician friend).

In JAGS/BUGS we can define the prior distribution on the precision of a normal distribution using a gamma distribution, which also happens to be a conjugate prior for the normal distribution parameterized by precision. We want to be able to specify a gamma prior on the normal distribution using our guess of the mean SD of the normal distribution and the SD of the SD of the normal distribution. In order to do this we need to find the prior distribution that corresponds to the gamma distribution but is the conjugate prior for a normal distribution parameterized by SD.

I found three mentions of this distribution where it is either called the inverted half gamma (Fink, 1997) or the inverted gamma-1 (Adjemian, 2010; LaValle, 1970).

The inverted gamma-1 distribution has two parameters $\nu$ and $s$ which corresponds to $2 \cdot shape$ and $2 \cdot rate$ of the gamma distribution respectively. The mean and SD of the inverted gamma-1 is:

$$ \mu = \sqrt{\frac{s}{2}}\frac{\Gamma(\frac{\nu-1}{2})}{\Gamma(\frac{\nu}{2})} \space \text{ and } \space \sigma^2 = \frac{s}{\nu - 2} - \mu^2$$

It doesn't seem to exist a closed solution that allow us to get $\nu$ and $s$ if we specify $\mu$ and $\sigma$. Adjemian (2010) recommends a numerical approach and fortunately a matlab script that does this is available from the open source platform Dynare. The following is an R translation of that script:

# Copyright (C) 2003-2008 Dynare Team, modified 2012 by Rasmus Bååth
#
# This file is modified R version of an original Matlab file that is part of Dynare.
#
# Dynare is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# Dynare is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.

inverse_gamma_specification <- function(mu, sigma) {
  sigma2 = sigma^2
  mu2 = mu^2
  if(sigma^2 < Inf) {
    nu = sqrt(2*(2+mu^2/sigma^2))
    nu2 = 2*nu
    nu1 = 2
    err = 2*mu^2*gamma(nu/2)^2-(sigma^2+mu^2)*(nu-2)*gamma((nu-1)/2)^2
    while(abs(nu2-nu1) > 1e-12) {
      if(err > 0) {
        nu1 = nu
        if(nu < nu2) {
          nu = nu2
        } else {
        nu = 2*nu
        nu2 = nu
        }
      } else {
        nu2 = nu
      }
      nu =  (nu1+nu2)/2
      err = 2*mu^2*gamma(nu/2)^2-(sigma^2+mu^2)*(nu-2)*gamma((nu-1)/2)^2
    }
    s = (sigma^2+mu^2)*(nu-2)
  } else {
    nu = 2
    s = 2*mu^2/pi
  }
  c(nu=nu, s=s)
}   

The R/JAGS script below shows how we can now specify our gamma prior on the precision of a normal distribution.

library(rjags)

model_string <- "model{
y ~ dnorm(0, tau)
sigma <- 1/sqrt(tau)
tau ~ dgamma(shape, rate)
}"

# Here we specify the mean and sd of sigma and get the corresponding
# parameters for the gamma distribution.
mu_sigma <- 100
sd_sigma <- 50
params <- inverse_gamma_specification(mu_sigma, sd_sigma)
shape <- params["nu"] / 2
rate <- params["s"] / 2

data.list <- list(y=NA, shape = shape, rate = rate)    
model <- jags.model(textConnection(model_string), 
    data=data.list, n.chains=4, n.adapt=1000)
update(model, 10000)
samples <- as.matrix(coda.samples(
    model, variable.names=c("y", "tau", "sigma"), n.iter=10000))

And now we can check if the sample posteriors (which should mimic the priors as we gave no data to the model, that is, y = NA) are as we specified.

mean(samples[, "sigma"])
## 99.87198
sd(samples[, "sigma"])
## 49.37357

par(mfcol=c(3,1), mar=c(2,2,2,2))
plot(density(samples[, "tau"]), main="tau")
plot(density(samples[, "sigma"]), main="sigma")
plot(density(samples[, "y"]), main="y")

Posterior plot

This seems to be correct. Any objections or comments to this method of specifying a prior is much appreciated!

Edit:

To calculate the shape and rate of the gamma prior as we have done above is not the same as directly using a gamma prior on the SD on a normal distribution. This is illustrated by the R script below.

# Generating random precision values and converting to 
# SD using the shape and rate values calculated above
rand_precision <- rgamma(999999, shape=shape, rate=rate)
rand_sd <- 1/sqrt(rand_prec)

# Specifying the mean and sd of the gamma distribution directly using the
# mu and sigma specified before and generating random SD values.
shape2 <- mu^2/sigma^2
rate2 <- mu/sigma^2
rand_sd2 <- rgamma(999999, shape2, rate2)

The two distributions now has the same mean and SD.

mean(rand_sd)
## 99.96195
mean(rand_sd2)
## 99.95316
sd(rand_sd)
## 50.21289
sd(rand_sd2)
## 50.01591

But they are not the same distribution.

plot(density(rand_sd[rand_sd < 400]), col="blue", lwd=4, xlim=c(0, 400))
lines(density(rand_sd2[rand_sd2 < 400]), col="red", lwd=4, xlim=c(0, 400))

enter image description here

From what I've read it seems to be more usual to put a gamma prior on the precision than a gamma prior on the SD. But I don't know what the argument would be for preferring the former over the latter.

References

Fink, D. (1997). A compendium of conjugate priors. http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.157.5540&rep=rep1&type=pdf

Adjemian, S. (2010). Prior Distributions in Dynare. http://www.dynare.org/stepan/dynare/text/DynareDistributions.pdf

LaValle, I.H. (1970). An introduction to probability, decision, and inference. Holt, Rinehart and Winston New York.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ @Procrastinator Yeah, I know that I'm just using rjags as a glorified random number generator. I though it would be useful to show as starting point for actually building a model in JAGS, and it is what I would have wanted to read when I first asked this question :) I've also added a comparison between your way of specifying the SD and my way in the end of my answer. Does it make sense? $\endgroup$ – Rasmus Bååth Nov 8 '12 at 13:39
  • $\begingroup$ I want to use the gamma on the precision of the normal distribution, as it is conjugate ,as you say, and often used in practice. But I don't want to specify the gamma prior on precision by it's shape and rate, as I feel it is non intuitive. I don't want to specify the gamma prior on precision with the mean and SD (of the precision). I want to specify the gamma using the mean of sd the normal distribution and sd of the sd of the normal distribution as it is easier. If the normal distribution was parameterized by sd and the "inverted half gamma" existed in JAGS I could have used them directly. $\endgroup$ – Rasmus Bååth Nov 8 '12 at 14:23
  • 1
    $\begingroup$ Gamma prior on the precision is what I want to use! My question (and answer) was about how to specify that gamma so that the resulting prior on the precision corresponded to a prior on the SD (1/sqrt(precision)) with a mean SD and SD of the SD that I specified my self. It's a shame that we could not meet IRL, communicating is so much easier face to face :) $\endgroup$ – Rasmus Bååth Nov 8 '12 at 14:53
  • 2
    $\begingroup$ I did this exact calculation a few years ago. If you're willing to credit me as someone-who-knows-what-he/she-is-doing, be reassured that your solution is correct. $\endgroup$ – Cyan Nov 8 '12 at 23:42
  • $\begingroup$ @Cyan Great to hear :) I was a bit afraid that there was something wrong with my question since it didn't get any up votes. What's your opinion on this reparameterization (I believe it could be quite useful myself)? Did you end up using it in practice or are there other alternatives? If you want to formulate this as an answer it would be possible to give you the bounty as it would go to waste otherwise... $\endgroup$ – Rasmus Bååth Nov 9 '12 at 7:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.