# Joint posterior distribution of differences

I have data $$x_1,...,x_n$$, $$y_1,...,y_m$$ and $$z_1,...,z_p$$ where $$x_1,...,x_n\sim N(\mu_x,\sigma^2_x)$$ and $$y_1,...,y_m\sim N(\mu_y,\sigma^2_y)$$ and $$z_1,...,z_p\sim N(\mu_z,\sigma^2_z)$$

Now let's assume I want to take a Bayesian approach and place the following priors: $$p(\mu_x,\sigma^2_x)\propto (\sigma^2_x)^{-1}$$, $$p(\mu_y,\sigma^2_y)\propto (\sigma^2_y)^{-1}$$, and $$p(\mu_z,\sigma^2_z)\propto (\sigma^2_z)^{-1}$$. Given these priors, I know what the posterior distribution is, but more importantly, I know that the conditional marginal distributions are

$$\mu_x|\sigma^2_x,x_1,...,x_n\sim N(\bar{x},\sigma^2_x/n)$$ and $$\mu_y|\sigma^2_y, y_1,...,y_m\sim N(\bar{y},\sigma^2_y/m)$$ and $$\mu_z|\sigma^2_z,z_1,...,z_n\sim N(\bar{z},\sigma^2_z/p)$$

where $$\bar x$$ is the average of the $$x$$'s. Similarly, for the case of the $$y$$'s and $$z$$'s.

I am interested in deriving the joint distribution of $$\mu_x-\mu_y$$, $$\mu_y-\mu_z$$, $$\mu_y-\mu_z$$ . Does the following approach make sense?

We know that the conditional posterior distribution of $$\mu_x$$ is $$p(\mu_x|\sigma^2_x,x_1,...,x_n) = N\left(\bar x, \frac{\sigma^2_x}{n}\right)$$

and similarly for $$p(\mu_y|\sigma^2_y,y_1,...,y_m)$$ and $$p(\mu_z|\sigma^2_z,z_1,...,z_p)$$.

Now, define \begin{align} \Delta := \begin{bmatrix} \delta_{xy}\\\\ \delta_{xz}\\\\ \delta_{yz} \end{bmatrix} = \begin{bmatrix} \mu_x-\mu_y\\\\ \mu_x - \mu_z\\\\ \mu_y - \mu_z \end{bmatrix} = \begin{bmatrix} 1 & -1 & 0\\\\ 1 & 0 & -1\\\\ 0 & 1 & -1 \end{bmatrix}\begin{bmatrix} \mu_x\\\\ \mu_y\\\\ \mu_z \end{bmatrix} =: A\mu, \end{align} then we have that $$\Delta|\sigma^2_x,\sigma^2_y,\sigma^2_z, x_1,...,x_n, y_1,...,y_m, z_1,...,z_p \sim N_3\left(A \mathbb{E}(\mu), A\mathbb{C}\text{ov}(\mu)A^T\right)$$ where \begin{align*} \mathbb{E}(\mu) = \begin{bmatrix} \bar{x}\\\\ \bar{y}\\\\ \bar{z} \end{bmatrix} \text{ and } \mathbb{C}\text{ov}(\mu) = \begin{bmatrix} \sigma^2_x/n & 0 & 0 \\\\ 0 & \sigma^2_y/m & 0 \\\\ 0 & 0 & \sigma^2_z/p \end{bmatrix}. \end{align*}

So my ultimate questions are the following:

1. Would the distribution above for $$\Delta$$ be the correct joint conditional posterior distribution of the differences?

2. And if so, would the appropriate strategy to obtain posterior samples from $$\Delta$$ be to first samples $$\sigma^2_x, \sigma^2_y,$$ and $$\sigma^2_z$$ from their joint posterior distribution?

3. And if so for number 2) does the joint posterior distribution of $$\sigma^2_x, \sigma^2_y,$$ and $$\sigma^2_z$$ have a closed form?

1. The conditional analysis is correct (assuming posterior conditional independence between $$\mu_x,\mu_y,\mu_z$$), as already indicated in the answer to your earlier question!
2. Simulating $$σ_x$$, $$σ_y$$, and $$σ_z$$ from their joint (marginal) posterior distribution, then $$\mu_x,\mu_y,\mu_z$$ from the corresponding conditional posterior is correct. However, if simulating from the (marginal) posterior proves complex, a Gibbs sampling approach may be preferable. But...
3. ...since $$\int_{-\infty}^{+\infty} \exp\left\{-\left[n(\bar x_n-\mu_x)^2+\sum_{i=1}^n(x_i-\bar x_n)^2\right]\big/2\sigma^2_x\right\}\,\frac{\text d\mu_x}{\sigma_x^{n+2}}\\\propto\exp\left\{-\sum_{i=1}^n(x_i-\bar x_n)^2\big/2\sigma^2_x\right\}\sigma_x^{-n-1}$$ and the equivalent closed forms for $$σ^2_y$$, and $$σ^2_z$$, the marginal posterior on $$(\sigma^2_x,\sigma^2_y,\sigma^2_z)$$ is available.
• Thanks for confirming 1) and 2). For 3) I understand that the marginal posteriors for the $\sigma^2$s are available in closed form (inverse-gammas), but are you implying that that joint (marginal) of the $\sigma^2$'s is not available in closed form (I was hoping that it reduced some how to, say, and inverse wishart)? Also, is there any reason or probability rule I am breaking by sampling from the marginals of each $\sigma^2$ and then plugging those into the joint distribution of $\Delta$ to obtain a posterior sample of $\Delta$ from the corresponding Normal distribution? Nov 19, 2022 at 14:55
• Sorry for the multiple questions, but if I interpret your answer correctly, that is to say I have $p(\Delta, \sigma^2_x, \sigma^2_y,\sigma^2_z|data) = p(\Delta| \sigma^2_x, \sigma^2_y,\sigma^2_z,data) p( \sigma^2_x, \sigma^2_y,\sigma^2_z|data)$ and here $p(\sigma^2_x, \sigma^2_y,\sigma^2_z|data) = p( \sigma^2_x|data)p( \sigma^2_y|data)p( \sigma^2_z|data)$ by independence? Nov 19, 2022 at 17:10
• Sorry, one more question. As I was going to code this up, I realize I don't recognize the kernel of the last expression in your 3). At first I had assumed it was inverse-gamma, however inside the exponential you have $\sigma^2_x$ and outside the exponential you have $\sigma_x^{-n-1}$. In order to get it to match upI assume a bit of reworking is needed. I want to confirm that the marginal of $\sigma^2_x$ is Inverse-gamma($n/2 - 1/2, \Sigma(x_i-\bar{x})^2/2$) Nov 20, 2022 at 16:35
• I get it now! I was forgetting about the original posterior distribution of $(\mu,\sigma^2)$ and was only thinking in terms of $\Delta$. Nov 23, 2022 at 14:40