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

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?

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.

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.

$\begingroup$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.$\endgroup$
– Mr. ZenMay 5 '18 at 16:38

$\begingroup$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.$\endgroup$
– BenMay 5 '18 at 23:55

$\begingroup$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?$\endgroup$
– Mr. ZenMay 7 '18 at 15:48