I'm looking for a resource that explains how to do update the coefficients for a Bayesian gamma regression using Gibbs sampling. Specifically, if

$y_i \sim Gamma(\alpha,\beta_i)$ and my data likelihood is $\prod_{i=1}^NP(y_i | \alpha,\beta_i)$

$log(\mu_i) = x_i \gamma_1 + \gamma_0$

$\alpha = \mu_i^2 / \phi$

$\beta_i = \mu_i/\phi$

assuming the shape parameter $\alpha$ is a known constant, I need to sample $\gamma_1$ and $\gamma_0$ from the full conditionals. For example, $P(\gamma_1 | \gamma_0,Y,X,\alpha) \propto P(y_i | \alpha, \alpha /exp(x_i \gamma_1 + \gamma_0))P(\gamma_1|\mu_0,\sigma^2_0)$.

However, the first term is gamma density and second is normal and are thus not conjugate. What is the best approach to sampling from this conditional? I can think of three options but I'm not sure what's best in terms of convergence time.

1) Metropolis proposals for $\gamma_1$. Proposal width can be problematic.

2) Adaptive rejection sampling. But I'm not sure about concavity.

3) Normal approximation for gamma. Can't seem to figure this out when variance and mean have coupled terms.

Note that I am aware I can just do this in STAN (or equivalent). However, I would still like to learn more about the theory behind this problem (i.e. not a application to real data, just want to know how it works). In fact, the reason I'm asking this question is because I can't seem to find any resources on how to do this besides "use STAN".


Your Answer

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

Browse other questions tagged or ask your own question.