I have been reading a couple related papers using Bayesian inference in hierarchical models1,2,3 but am struggling to bridge the gap in one aspect of the papers. I think the struggle is in relation to the posterior predictive distribution. The model is described as

$$log(y_{i,t}) \sim \mathcal{N}(\beta_{0,i} + \beta_{1,i}a_{i,t} + \eta_t, \sigma_y^2)$$

$$\eta_t \sim \mathcal{N}(\beta_2 x_t, \sigma_{\eta}^2)$$

$$x_t \sim \mathcal{N}(\mu_x, \sigma_x^2)$$

In this case, $y_{i,t}$, $a_{i,t}$, $x_t$ are measured but the goal will be to predict new values of $x_t$ (climate) for which we have measures of $y_{i,t}$ and $a_{i,t}$. They state that the posterior predictive distribution can be sampled from

$$x_t^{(j)} \sim \mathcal{N}\left(\frac{\sigma_{\eta}^{2(j)}\mu_x^{(j)} + \sigma_x^{2(j)}\beta_2^{(j)}\eta_t^{(j)}}{\sigma_{\eta}^{2(j)} + \sigma_x^{2(j)}\beta_2^{2(j)}}, \left[\frac{1}{\sigma_x^{2(j)}} + \frac{\beta_2^{2(j)}}{\sigma_{\eta}^{2(j)}} \right] \right) $$

where $(j)$ represents the $j^{th}$ MCMC sample. I know that the posterior predictive distribution is defined as

$$p(\tilde{x} \mid x) = \int_\theta p(\tilde{x} \mid \theta)p(\theta \mid x)d\theta$$

However, I am unable to get from the model description to the posterior using this equation. Could anyone walk me through the probability/integration steps necessary to come up with this specific posterior predictive distribution?

1. Schofield et al. 2016

2. Steinschneider et al. 2017

3. Schofield and Barker 2017

  • $\begingroup$ Let me know if this would be better posted on math.stackexchange.com $\endgroup$ – djhocking May 6 '20 at 18:23

It appears that the stated distribution is for $x_t^{(j)} | \eta_t^{(j)}$ and the random variable $y_{i,t}$ is being ignored for now. It also appears that the authors are being a bit loose in their notation for the normal distribution, using the variance parameter in some statements and the precision parameter in others. (I will parameterise with the variance unless otherwise stated.) To obtain the conditional density we take the joint density kernel and "complete the square" to simplify. Taking proportionality with respect to the $x$ variable gives:

$$\begin{aligned} p(x| \eta) &\equiv p(x_t^{(j)} = x| \eta_t^{(j)} = \eta) \\[12pt] &\overset{x}{\propto} p(x_t^{(j)} = x, \eta_t^{(j)} = \eta) \\[12pt] &= p(\eta_t^{(j)} = \eta | x_t^{(j)} = x) \cdot p(x_t^{(j)} = x) \\[12pt] &= \text{N}(\eta | \beta_2 x, \sigma_\eta^2) \cdot \text{N}(x | \mu_x, \sigma_x^2) \\[6pt] &\overset{x}{\propto} \exp \bigg( - \frac{1}{2 \sigma_\eta^2} (\eta - \beta_2 x)^2 \bigg) \cdot \exp \bigg( - \frac{1}{2 \sigma_x^2} (x - \mu_x)^2 \bigg) \\[6pt] &= \exp \bigg( - \frac{1}{2} \bigg[ \frac{1}{\sigma_\eta^2} (\eta - \beta_2 x)^2 + \frac{1}{\sigma_x^2} (x - \mu_x)^2 \bigg] \bigg) \\[6pt] &= \exp \bigg( - \frac{1}{2} \bigg[ \frac{1}{\sigma_\eta^2} (\eta^2 - 2 \eta \beta_2 x + \beta_2^2 x^2) + \frac{1}{\sigma_x^2} (x^2 - 2 \mu_x x + \mu_x^2) \bigg] \bigg) \\[6pt] &= \exp \bigg( - \frac{1}{2} \bigg[ \Big( \frac{1}{\sigma_x^2} + \frac{\beta_2^2}{\sigma_\eta^2} \Big) x^2 -2 \Big( \frac{\mu_x}{\sigma_x^2} + \frac{\eta \beta_2}{\sigma_\eta^2} \Big) x + \Big( \frac{\eta^2}{\sigma_\eta^2} + \frac{\mu_x^2}{\sigma_x^2} \Big) \bigg] \bigg) \\[6pt] &= \exp \bigg( - \frac{1}{2} \Big( \frac{1}{\sigma_x^2} + \frac{\beta_2^2}{\sigma_\eta^2} \Big) \bigg[ x^2 -2 \Big( \frac{\eta \sigma_x^2 \beta_2^2 + \mu_x \sigma_\eta^2}{\sigma_x^2 \beta_2^2 + \sigma_\eta^2} \Big) x + \text{const} \bigg] \bigg) \\[6pt] &\overset{x}{\propto} \exp \bigg( - \frac{1}{2} \Big( \frac{1}{\sigma_x^2} + \frac{\beta_2^2}{\sigma_\eta^2} \Big) \bigg( x - \frac{\eta \sigma_x^2 \beta_2^2 + \mu_x \sigma_\eta^2}{\sigma_x^2 \beta_2^2 + \sigma_\eta^2} \bigg)^2 \bigg) \\[6pt] &\overset{x}{\propto} \text{N}\bigg( x \Bigg| \text{Mean} = \frac{\eta \sigma_x^2 \beta_2^2 + \mu_x \sigma_\eta^2}{\sigma_x^2 \beta_2^2 + \sigma_\eta^2}, \text{Precision} = \frac{1}{\sigma_x^2} + \frac{\beta_2^2}{\sigma_\eta^2} \bigg). \\[6pt] \end{aligned}$$

Thus, the conditional distribution is:

$$x_t^{(j)} | \eta_t^{(j)} \sim \text{N}\bigg(\text{Mean} = \frac{\eta_t^{(j)} \sigma_x^2 \beta_2^2 + \mu_x \sigma_\eta^2}{\sigma_x^2 \beta_2^2 + \sigma_\eta^2}, \text{Precision} = \frac{1}{\sigma_x^2} + \frac{\beta_2^2}{\sigma_\eta^2} \bigg).$$

  • $\begingroup$ That's good to know that they changed to using the precision in this case. They were using JAGS for implementing their Bayesian model which is parameterized in terms of the precision, which is probably why they switched. I also believe the third paper (a response) criticizes the use of the conditional distribution that ignores y. Thanks for the detailed answer. $\endgroup$ – djhocking May 10 '20 at 1:22
  • $\begingroup$ Glad my answer helped. In Bayesian statistics it is usually more convenient to parameterise the normal distribution in terms of the precision, rather than the variance (for reasons such as the above). $\endgroup$ – Ben May 10 '20 at 2:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.