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

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

and

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

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{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}$$

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. There are two publications that deal with a similar problem (however, without any guidance on the integration process). There, the marginal prior includes the modified Bessel function of the second kind and the posterior moments of the $GIG$ distribution. That's why I think that this problem might not be trivial and I'm thankful for help.

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{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}

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? – Mr. Zen May 4 '18 at 3:46
• Yeah, that doesn't really help you. – Ben May 4 '18 at 3:50
• 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. – Mr. Zen May 5 '18 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. – Ben May 5 '18 at 23:55
• 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? – Mr. Zen May 7 '18 at 15:48