Posterior predictive distribution of the difference

Question

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.

Answer 1

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)$$