Bayesian Hierarchical Model - Exact Conjugate Solution? I was hoping to get some help. In understand how to compute an exact numerical solution (http://www.cs.berkeley.edu/~jordan/courses/260-spring10/lectures/lecture5.pdf) for the following Bayesian model:
$$ \tau \sim Ga(\alpha, \beta)$$
$$\mu \sim N(m ,p)$$
$$Y_i \sim N(\mu, \tau)$$
where the data is $Y_i$.
I am performing a meta-analysis where my data are $Y_i$ from each study, and $\tau_i$ from each study. The problem is described by a hierarchical model:
$$ \tau \sim Ga(\alpha, \beta)$$
$$\mu \sim N(m ,p)$$
$$\theta_i \sim N(\mu, \tau)$$
$$Y_i \sim N(\theta_i, \tau_i)$$
Is there any way to solve this model analytically? Any help/suggestions would be greatly appreciated.
 A: I'm assuming that by "exact numerical solution" and "analytically" you mean that you want $$p(\mu,\theta_i,\tau|Y_1,\ldots,Y_n, \tau_1,\ldots,\tau_n)$$ to be a known distribution. If this is what you want, then NO there is no way to solve for this posterior analytically. 
Let me use $\sigma^2$ instead of your $\tau$, so that we have $$\theta_i \stackrel{iid}{\sim} N(\mu,\sigma^2)$$. 
Some suggestions that would get you closer to an analytic solution are 


*

*Let $Y_i\stackrel{iid}{\sim} N(\theta_i,\sigma^2\tau_i)$.

*Let $\mu \sim N(m,\sigma^2p)$ 

*Let $\sigma^2 \sim IG(\alpha,\beta)$. 


Now you should be able to integrate out $\theta_1,\ldots,\theta_n,\mu$ to obtain $Y\sim N(m, \sigma^2 S)$ where $Y=(Y_1,\ldots,Y_n)$and $S$ is a known matrix. Then, you should be able to get $\mu|Y,\tau_1,\ldots,\tau_n,\sigma^2$ followed by $\theta_1,\ldots,\theta_n|Y,\tau_1,\ldots,\tau_n,\mu,\sigma^2$. Although this wouldn't get you to an analytic solution, you can at least use Monte Carlo (rather than Markov chain Monte Carlo) to obtain draws from the full posterior. 
A: Unfortunately the algebraic structure of the normal likelihood does not allow us to separate the $\tau_i$ and $\tau$ in a way that allows for exact conjugate relationships to be used here. The exponential family conjugate relationships are a direct consequence of the sum/product properties of exponentials.. to see the problem look at the log likelihood of the data:
$$
\text{LL}(\text{data}) = \text{constant} + \frac{1}{2}\sum_i \log(\tau_i) +  \frac{1}{2}\sum_i \tau_i (Y_i - \theta_i)^2.
$$
There is no way to combine terms involving $\theta_i$ with the prior for $\theta_i$,
$$
\log(p(\theta_i)) = \text{constant} + \frac{1}{2} \log(\tau) + \frac{1}{2} \tau (\theta_i - \mu)^2.
$$
To combine these two (and the other distributions involved naturally) and get a posterior distribution as a function of the parameters, we would normally combine the squared sum at the end of each of these through completing the square (see section 3 here). This is not possible here because each of the $\tau_i$ infront of each $(Y_i - \theta_i)^2$ are not the same. By not being able to complete the square we are not able to put  the parameters in a normal distribution form, so there is by definition no conjugacy. This is without even investigating the additional levels of depth necessitated by the prior (it's really just one Normal-Gamma prior) on $\tau$ and $\mu$.
Further, this non-conjugacy result would hold if you use the convolution of normals mentioned by SeanEaster due to the way the variances are added, making them impossible to wrangle into a Normal-Gamma posterior.

An alternative approach to working with a similar model specification is given by jaradniemi
