The hierarchical model specified below is quite "standard" and easy to implement in for example JAGS/BUGS. It has an hierarchical gamma prior on the subject's precisions ($\tau_j$) which in turn has gamma priors on its mean ($m$) and SD ($d$). After having fitted this model I have access to estimates of the mean group precision ($m$) and SD of the group precision ($d$) but I rather make inferences on the natural scale. That is, the mean group SD and the SD of the group SD.

How can I make inferences on the mean group SD and the SD of the group SD from such a model? Is it at all possible?

Here is a Kruschke style diagram of the model taken from the book Doing Bayesian Data Analysis

enter image description here

And here is the corresponding BUGS code:

model {
  for( i in 1 : Ndata ) {
    y[i]  ̃ dnorm( mu[ subj[i] ] , tau[ subj[i] ] )
  for ( j in 1 : Nsubj ) {
    mu[j]  ̃ dnorm( muG , tauG )
    tau[j]  ̃ dgamma( sG , rG )
  muG  ̃ dnorm( 2.3 , 0.1 )
  tauG  ̃ dgamma( 1 , .5 )
  sG <- pow(m,2) / pow(d,2)
  rG <- m / pow(d,2)
  m  ̃ dgamma( 1 , .25 )
  d  ̃ dgamma( 1 , .5 )

One solution would perhaps be to not put the hyperpriors on the precision ($\tau_j$) but rather on the SD ($\sigma_j$) though I'm not at all sure what priors I should use then.


1 Answer 1


Given that the $\tau_j$ are i.i.d. Gamma with parameters $\alpha, \lambda$ (where $f(x) \propto \lambda^{\alpha} x^{\alpha-1}\exp\{-\lambda x\}$, to define the parameterization of the Gamma), then $\sigma_j = 1/\sqrt{\tau_j}$ has a Stacy distribution, which is a special case of the Amoroso distribution. This has a closed form expression for the mean and standard deviation, assuming the parameters are such that the mean and standard deviation exist, in terms of $\alpha$ and $\lambda$, and hence in terms of your m and d.

Define $\theta = 1/\sqrt{\lambda}$. Then:

$\mathbb{E}\sigma_j = \alpha + \theta {\Gamma(\alpha-{1\over 2}) \over \Gamma(\alpha)} $


$\text{Sd}\space \sigma_j = \theta \left[{\Gamma(\alpha-1) \over \Gamma(\alpha)}-\left({\Gamma(\alpha-{1\over 2}) \over \Gamma(\alpha)}\right)^2 \right]^{1\over2}$

where $\alpha > 1/2$ for the mean to exist and $\alpha > 1$ for the variance to exist. (If you read the article, note the $\beta$ parameter of the distribution $= -2$, as $\tau = \sigma^{-2}$.)

One method for generating MCMC samples of the mean and std. deviation of the $\sigma_j$ would be to take your samples of m and d (or, more simply, rG (= $\lambda$) and sG (= $\alpha$)), and transform them as above. However, one suspects that you may indeed have samples where $\alpha \leq 1$ or even $\alpha \leq 1/2$. No doubt many solutions to this exist, but two that come to mind are to 1) throw out all the samples (in their entirety) where this happens, which corresponds to altering your prior on m and d so that they are confined to a region of $\mathbb{R}^{+2}$ where the mean and variance of the $\tau$ exist, or 2) put your prior on $\sigma$ instead of $\tau$, altering the hyperpriors to something appropriate. Without knowing how strongly you're attached to the particular shape of the priors on m and d, though, it's hard to give advice about what "something appropriate" might be. Maybe just put a Gamma distribution on the $\sigma_j$ instead of the $\tau_j$ and continue as you did... hopefully the choice of prior doesn't make much difference to your posterior.

If it doesn't happen all that often, and you like those priors as they are, I'd probably go with option 1. If it happens a lot, though, then those priors in combination with your data are perhaps telling you something about the upper tail of the posterior distribution of the $\sigma_j$.

To summarize why this works, this is equivalent to:

1) Put a Stacy prior on the $\sigma_i$, parameterized by rG and sG,

2) Put hyperpriors on rG and sG as they currently are,

3) Instead of bothering to generate Stacy-distributed $\sigma_i$ and transforming them to $\tau_i$, just generate the $\tau_i$ directly using the relationship between the Stacy prior on $\sigma$ and the equivalent Gamma prior on $\tau$,

4) Calculate samples of the first two moments of the Stacy prior from the sampled rG and sG.

  • $\begingroup$ Thank you! I've been looking for this answer for so long :'-) Really, thanks! You are correct in that I get some (actually a lot of) samples where α < 1. While I try to figure out why that is the case I would be happy to put a prior on σ meanwhile. The reason why I haven't done so yet is because I haven't found any example of someone putting informative on σ, it is always just gamma on τ. What I'm looking for is a reasonable distribution for σ that can be parameterized with scale and location and that is not too "unconventional". Would you say the gamma distribution could work in this case? $\endgroup$ Nov 28, 2012 at 21:49
  • 1
    $\begingroup$ I know you've seen this link before, but stats.stackexchange.com/questions/6493/… can help; many of the mentioned distributions can be parameterized by a mean and variance which are transformed to the "usual" parameters (as you did with the Gamma above). As to why you get many $\alpha < 1$, it's because the tails of your prior on $\sigma$, when you integrate out m and d, are too fat for $\sigma$ to have a finite variance. You might also be able to fix it by bounding $d$ above and below (JAGS lets you do this.) $\endgroup$
    – jbowman
    Nov 28, 2012 at 21:59
  • $\begingroup$ Yes, I looked at that page many times :) Again! $\endgroup$ Nov 28, 2012 at 22: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.