# Deriving posteriors in a hierarchical model

So I'm doing a problem with the following set up:

$$X | \mu_1 \sim N(\mu_1, \sigma_1^2)$$ $$\mu_1 | \mu_2 \sim N(\mu_2, \sigma_2^2)$$ $$...$$ $$\mu_{k-1} | \mu_k \sim N(\mu_k, \sigma_k^2)$$

where $$\sigma_i^2$$ are all known. The problem is to derive the posteriors $$\pi(\mu_i|X, \mu_k)$$.

My attempt has been to first derive the joint density as the product of conditionals

$$f(X,\mu_1,\dots,\mu_{k-1}) = f(X|\mu_1)\prod_{i=1}^{k-1} f(\mu_i |\mu_{i+1})$$

and then from here to integrate out $$\mu_j \neq \mu_i$$. This leaves us with the joint distribution $$f(X,\mu_i,\mu_k)$$ and so proportionality can be used to identify the posterior distribution $$f(\mu_i |X,\mu_k)$$.

However, the integral is not simple, and I have not made any progress. I'm wondering if there's a simpler way to go about this? I know the answer, which is a given normal distribution, but getting there is more difficult.

$$\mu_{k-1} = \mu_k + \epsilon_k$$ $$\mu_{k-2} = \mu_k + \epsilon_k + \epsilon _{k-1}$$ $$...$$ $$\mu_{i} = \mu_k + \epsilon_k + \epsilon_{k-1} + ...+\epsilon_{i+1}$$ $$...$$ $$X = \mu_k + \epsilon_k + \epsilon_{k-1} + ...+\epsilon_i + ...+ \epsilon_1$$
where $$\epsilon_i \sim N(0,\sigma_i^2)$$
Because $$X$$ and $$\mu_i$$ are linear combinations of normal distributed random variables, so they are bivariate normal distributed. It is easy to get their joint distribution, then get the conditional distribution.