0
$\begingroup$

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.

$\endgroup$
3
  • 1
    $\begingroup$ The second strategy works just fine. For justification google search “law of the unconscious statistician” $\endgroup$
    – Taylor
    Dec 1, 2022 at 18:42
  • $\begingroup$ @Taylor, I updated my question. Do you think the last strategy is correct? $\endgroup$
    – John Smith
    Dec 1, 2022 at 18:59
  • 1
    $\begingroup$ 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. $\endgroup$
    – Taylor
    Dec 1, 2022 at 21:28

1 Answer 1

1
$\begingroup$

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

$\endgroup$
7
  • 1
    $\begingroup$ 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? $\endgroup$
    – John Smith
    Dec 1, 2022 at 19:57
  • $\begingroup$ 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$. $\endgroup$
    – John Smith
    Dec 1, 2022 at 20:02
  • $\begingroup$ 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$ $\endgroup$
    – John Smith
    Dec 1, 2022 at 21:00
  • 1
    $\begingroup$ I think I get it then. If you can write one of the components of the random vector as as a linear combination of the others, then its a degenerate normal, correct? $\endgroup$
    – John Smith
    Dec 2, 2022 at 22:40
  • 1
    $\begingroup$ Actually, I just realized that your comment is more of a caution on how to sample from this multivariate distribution, but not that it is problematic (I read this the wrong way!). So I'm fine to use a degenerate multivariate normal as long as aI can sample from it correctly. And in this case I am using R's mvrnorm package which seems capable of doing so. $\endgroup$
    – John Smith
    Dec 2, 2022 at 23:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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