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.
 A: 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.
