# Marginal prior derivation in hierarchical Bayesian model

I am working on a model that is closely related to the normal gamma shrinkage prior setup discussed in Griffin & Brown (2010). Suppose we want to draw $$P$$ parameters $$\beta_p$$ with $$p=1,...,P$$. The prior setup is given by

$$\begin{equation} \beta_{p} \sim N(0, \frac{2}{\lambda}\tau^2_{p}) \end{equation}$$

and

$$\begin{equation} \begin{split} \tau^2_{p} &\sim G(\theta_1, \theta_2) \\ \lambda &\sim G(c_0,c_1) \end{split} \end{equation}$$

where $$G$$ denotes the gamma distribution. This is usually refered to as global/local shrinkage prior setup where $$\tau^2_{p}$$ is the local shrinkage parameter and $$\lambda$$ is the global shrinkage parameter.

I would like to derive the marginal prior for $$\beta$$, but do not really know, how to proceed. Can anyone point me in the right direction?

The conditional posterior distributions can be derived as

$$\begin{equation} \begin{split} g_1 &= P * \theta + c_0\\ d &= c_1 + \frac{\theta}{2\sum_{p=1}^P \tau^2_p}\\ \ \\ \pi(\lambda) &\sim G(g_1,d)\\ \pi(\tau^2_{p}~|~\cdot) &\sim GIG(\theta-0.5,~ \beta_{p}^2,~ \lambda \theta)\\ \end{split} \end{equation}$$

where $$GIG$$ denotes the Generalized Inverse Gaussian distribution.

I know I'm supposed to integrate out $$\lambda$$ and $$\tau_p^2$$, but I have no idea how to start.

References

Griffin, J. E., & Brown, P. J. (2010). Inference with normal-gamma prior distributions in regression problems. Bayesian Analysis, 5(1), 171-188.

Using the scaling rules for gamma random variables, we can write:

$$\tau_p^2 \sim \frac{\text{Chi-Sq}(2 \theta_1)}{2\theta_2} \quad \quad \quad \lambda \sim \frac{\text{Chi-Sq}(2 c_0)}{2 c_1}.$$

Hence, the ratio has a scaled F-distribution:

$$\phi_p \equiv \frac{\lambda}{\tau_p^2} \cdot \frac{\theta_1 c_1}{\theta_2 c_0} \sim \text{F}(2 c_0, 2 \theta_1).$$

The marginal prior for $\beta_p$ is given by the following integral:

\begin{equation} \begin{aligned} \pi(\beta_p) &= \int \limits_0^\infty \pi( \beta_p | \phi_p = r ) \pi( \phi_p = r ) dr \\[8pt] &= \int \limits_0^\infty \text{N} \Big( \beta_p \Big| 0, \frac{\theta_1 c_1}{\theta_2 c_0} \frac{2}{r} \Big) \text{F} ( r | 2 c_0, 2 \theta_1 ) dr \\[8pt] &\propto \int \limits_0^\infty r^{1/2} \exp \Big( - \frac{\theta_2 c_0}{\theta_1 c_1} \beta_p^2 \cdot r \Big) r^{c_0-1} \Big( 1+ \frac{c_0}{\theta_1} r \Big)^{-c_0-\theta_1} dr \\[8pt] &= \int \limits_0^\infty \exp \Big( - \frac{\theta_2 c_0}{\theta_1 c_1} \beta_p^2 \cdot r \Big) r^{c_0-1/2} \Big( 1+ \frac{c_0}{\theta_1} r \Big)^{-c_0-\theta_1} dr \\[8pt] &\propto \int \limits_0^\infty \exp \Big( - \frac{\theta_2}{\theta_1} \frac{2 c_0 + 1}{2 \theta_1 - 1} \beta_p^2 \cdot s \Big) s^{c_0-1/2} \Big( 1+ \frac{2 c_0 + 1}{2 \theta_1 - 1} s \Big)^{-c_0-\theta_1} ds \\[8pt] &\propto \int \limits_0^\infty \exp \Big( - \frac{\theta_2}{\theta_1} \frac{2 c_0 + 1}{2 \theta_1 - 1} \beta_p^2 \cdot s \Big) \text{F} (s | 2 c_0 + 1, 2 \theta_1 - 1) ds \\[8pt] &\propto \mathbb{E} \Bigg[ \exp \Big( - \frac{\theta_2}{\theta_1} \frac{2 c_0 + 1}{2 \theta_1 - 1} \beta_p^2 \cdot S \Big) \Bigg] \quad \quad \quad S \sim \text{F} (2 c_0 + 1, 2 \theta_1 - 1). \\[8pt] \end{aligned} \end{equation}

From this result we can see that the marginal density for $\beta_p$ is proportional to the expected value of a certain exponential-decay function for an F-distributed random variable. The integral required to compute this has no closed form solution, so you would need to estimate the density either by numerical integration or by simulating the expected value result.

• Wow, thank you very much for working this out! I suppose $\theta_1 = \theta_2$ and $c_0 = c_1$ does not really make things easier, does it? May 4, 2018 at 3:46
• Dear Ben, thanks again for your guidance in the process of working this out, it was very helpful. However, one question remains: If I understand it correctly, this implies that the marginal density can only be evaluated for $\theta_1 > 0.5$? Typical values for $\theta_1$ range between 0.1 and 0.5, therefore this is problematic. Any thoughts on that? Thank you. May 5, 2018 at 16:38
• That is not necessarily a problem - the last two steps above implicitly assume $\theta_1 > 0.5$, but if you don't have this then the previous integral results still hold.
• However, the previous results are divergent integrals in case $\theta_1<0.5$ as the exponential term will explode as the negative sign vanishes, right? Do you know of any trick to avoid the resulting overflow when integrating numerically? May 7, 2018 at 15:48