# Posterior predictive distribution of the difference

My question here is related to the following question I posted here: Joint posterior distribution of differences

With respect to that last question, what I want to discuss is how to appropriately sample from the posterior predictive distribution of the differences, i.e., the distribution $$p(x^*-y^*, x^*-z^*,y^*-z^*|x,y,z)$$. I assume to get this distribution I would need to calculate something like $$p(x^*-y^*, x^*-z^*,y^*-z^*|x,y,z)=\\ \int_\Delta\int_{\sigma^2_x}\int_{\sigma^2_y}\int_{\sigma^2_z}p(\Delta,\sigma^2_x,\sigma^2_y,\sigma^2_z| x, y, z)p(x^*-y^*, x^*-z^*,y^*-z^*|\Delta,\sigma^2_x,\sigma^2_y,\sigma^2_z x, y, z)d\Delta d\sigma^2_xd\sigma^2_yd\sigma^2_z$$

I know, as a general strategy in the case of a general posterior predictive distribution, say $$p(x^*|x) = \int_\Theta p(x^*|\theta,x)p(\theta|x)dx,$$ if you want to sample from $$p(x^*|x)$$, one strategy is to first sample a posterior draw $$\theta$$ and then to plug that $$\theta$$ into $$p(x^*|\theta,x)$$ and then sample an $$x^*$$ (i.e., a posterior predictive draw) from $$p(x^*|\theta,x)$$ which is just the same type of distribution as the likelihood.

Now returning to my question, I figured I could first sample $$\sigma^2_x$$, then sample $$\sigma^2_y$$, then sample $$\sigma^2_z$$, and then sample $$\Delta$$ (i.e., get posterior samples), and then to plug those posterior samples into $$p(x^*-y^*, x^*-z^*,y^*-z^*|\Delta,\sigma^2_x,\sigma^2_y,\sigma^2_z| x, y, z)$$ and to take a sample from that likelihood to get a posterior predictive draw. However, where I am stumped, is what the form of the likelihood $$p(x^*-y^*, x^*-z^*,y^*-z^*|\Delta,\sigma^2_x,\sigma^2_y,\sigma^2_z x, y, z)$$ actually is.

My alternative thought was to sample independently from $$p(x^*|x)$$, $$p(y^*|y)$$, and $$p(z^*|z)$$ (which I know how to do) and to then subtract those samples, i.e., $$x^*-y^*$$, $$x^*-z^*$$, and $$y^*-z^*$$, but my concern is aren't the those samples $$x^*-y^*$$, $$x^*-z^*$$, and $$y^*-z^*$$ independent of one another, while a sample of $$p(x^*-y^*, x^*-z^*,y^*-z^*|\Delta,\sigma^2_x,\sigma^2_y,\sigma^2_z x, y, z)$$ is not?

Or could I use the following strategy, first sample $$\sigma^2_x$$, then sample $$\sigma^2_y$$, then sample $$\sigma^2_z$$, and then sample $$\mu_1$$, $$\mu_2$$, and $$\mu_3$$ (i.e., get posterior samples), and then take a sample from $$\begin{pmatrix}x^*-y^*\\ x^*-z^*\\ y^*-z^*\end{pmatrix}\sim N_3\left(A\begin{pmatrix}\mu_1\\ \mu_2\\ \mu_3\end{pmatrix}, A\begin{pmatrix}\sigma^2_x & 0 &0\\ 0 & \sigma^2_y & 0\\ 0 & 0 & \sigma^2_z\end{pmatrix}A^T\right)$$

where $$\begin{pmatrix}\sigma^2_x & 0 &0\\ 0 & \sigma^2_y & 0\\ 0 & 0 & \sigma^2_z\end{pmatrix}$$ is the associated covariance matrix, and $$A = \begin{pmatrix}1& -1 &0\\ 1 & 0 & -1\\ 0 & 1 & -1\end{pmatrix}$$ is a matrix of contrasts.

• The second strategy works just fine. For justification google search “law of the unconscious statistician” Dec 1, 2022 at 18:42
• @Taylor, I updated my question. Do you think the last strategy is correct? Dec 1, 2022 at 18:59
• You didn't mention the last strategy or anything about normality the last time I read this. Looks like you have your question answered, anyway, so I won't go any further. Dec 1, 2022 at 21:28

The predictive distribution$$p(x^*-y^*, z^*-x^*,y^*-z^*|\Delta,\sigma^2_x,\sigma^2_y,\sigma^2_z x, y, z)$$is a Normal distribution as a linear transform of a Normal vector: $$\delta\equiv\left[\begin{matrix} x^*-y^*\\ z^*-x^*\\y^*-z^* \end{matrix}\right] =\underbrace{\left[\begin{matrix} 1 &-1 &0\\-1 & 0& 1\\0 &1 &-1 \end{matrix}\right]}_{\mathbf D} \times \left[\begin{matrix} x^*\\ y^*\\z^* \end{matrix}\right]$$ Hence $$\delta=\left[\begin{matrix} \delta_1\\\delta_2\\\delta_3 \end{matrix}\right]\sim\mathcal N\left( \left[\begin{matrix} \mu_x-\mu_y\\ \mu_z-\mu_x\\\mu_y-\mu_z \end{matrix}\right], \mathbf D\, \text{diag}(\sigma^2_x,\sigma^2_y,\sigma^2_z)\,\mathbf D^\top \right)\tag{1}$$ but this Normal is degenerate since $$\delta_1+\delta_2+\delta_3=0$$ with probability one. To generate from (1), one thus need (only) generate $$(\delta_1+\delta_2)$$ and derive $$\delta_3=-\delta_1-\delta_2$$.
Since $$\mathbf D\, \text{diag}(\sigma^2_x,\sigma^2_y,\sigma^2_z)\,\mathbf D^\top=\left[\begin{matrix} \sigma^2_x+\sigma^2_y &-\sigma^2_x &-\sigma^2_y\\ -\sigma^2_x &\sigma^2_x+\sigma^2_z &-\sigma^2_z\\ -\sigma^2_y &-\sigma^2_z &\sigma^2_x+\sigma^2_y \end{matrix}\right]$$ the marginal of the pair is $$\left[\begin{matrix} \delta_1\\ \delta_2\\ \end{matrix}\right]\sim\mathcal N\left( \left[\begin{matrix} \mu_x-\mu_y\\ \mu_z-\mu_x \end{matrix}\right],\left[\begin{matrix} \sigma^2_x+\sigma^2_y &-\sigma^2_x\\ -\sigma^2_x &\sigma^2_x+\sigma^2_z \end{matrix}\right]\right)$$
• Are the variances in your covariance matrix correct? Shouldn't you take $Cov([x^*, y^*,z^*]^T)$ which is a diagonal matrix of individual variances rather than sums of two variances? Dec 1, 2022 at 19:57
• Thanks for clarifying! Also, can you fix the degenerate issue by simply changing $\delta_2$ to $\delta_2 = \mu_x-\mu_z$? You currently have it as $\delta_2 = \mu_z-\mu_x$. Dec 1, 2022 at 20:02
• I understand in the case you have shown its degenerate since knowing 2 of the random variables tells you the third, but in the case I describe, i.e., $d_2=\mu_x-\mu_z$ I don't understand why this is the case. If I took $\delta_1 + \delta_2 + \delta_3$ it does not sum to 0, i.e., $\delta_1 + \delta_2 + \delta_3= 2\mu_x-2\mu_z$ Dec 1, 2022 at 21:00