# Distribution of posterior mean from different datasets

This question has originated from this question.

Suppose we have the following simple setup, for $i = 1, \dots, n$ $$y_i \mid \mu \sim N( \mu, 1) \text{ and } \mu \sim N(0,1).$$

Then due to the nice conjugate setup the posterior distribution for $\mu$ is $$\mu \mid \mathbf{y} \sim N\left(\dfrac{\sum_{i=1}^{n} y_i}{n+1}, \dfrac{1}{n+1} \right).$$

Now suppose I do the following $B = 100$ times.

1. Simulate $n = 1000$ from $y \mid \mu$.
2. Draw $r = 2000$ samples from the posterior $X_1, \dots, X_r$ (no need for MCMC).
3. Calculate the posterior mean and posterior standard deviation est_mu and sd_mu.

So the code for this would look something like.

r <- 2000
B <- 100
n <- 1000
est_mu <- numeric(length = B)
sd_mu <- numeric(length = B)
x <- matrix(0, nrow = r, ncol = B)

for(i in 1:B)
{
# Sampling new data from likelihood
y <- rnorm(n, mean = 1, sd = 1)

# Drawing from the posterior
x[,i] <- rnorm(r, mean = sum(y)/(n+1), sd = sqrt(1/(n+1)))

est_mu[i] <- mean(x[,i])
sd_mu[i] <- sd(x[,i])
}
mean(sd_mu)
#[1] 0.03159523
sd(est_mu)
#[1] 0.03150087


The value for mean(sd_mu) makes sense as it is roughly $\sqrt{1/1001}$, but why is sd(est_mu) also roughly the same?

Question: How do you show that sd(est_mu) $\approx$ mean(sd_mu)?

Originally I thought that this was impossible, because the posterior mean for each of $j = 1, \dots B$ draw of the likelihood is $$\hat{\mu}_j = \dfrac{1}{r} \sum_{t=1}^{r} X_{r} \, \, \text{ where } X_r \sim \mu \mid Y_j$$

and thus the variance should be $\sigma^2/r$ where $\sigma^2 = Var(\mu|Y) = 1/1001$. But this is clearly not the case as witnessed in the above code.

If however, I use the same data everytime, I get the expected result of the variance being $\sigma^2/r$. So clearly, the randomness in generating the data is getting transferred to the posterior mean estimates and this makes sense intuitively. I am just not able to show it by math.

In addition, the OP on the original comment said that once you have the posterior means from $B = 100$ independent samples, this a just a sample from the posterior distribution, that is, $\hat{\mu}_i \sim \mu|Y$, and I am not able to show this result.

So maybe a more general question is

Question: Why is $\hat{\mu}_j \sim \mu|Y$?

• Your R code fails because object 'r' not found. It is also unclear why you have mean = 5 – Henry May 16 '16 at 6:37
• Presumably by $\mu \sim N (0, 1)$ you are referring to a prior rather than a distribution on a higher level of a hierarchical set-up? Is your result just the result of the particular choice of true value of $\mu=5$ that you simulate with? It seems like you had a pretty bad prior... – Björn May 16 '16 at 7:52
• For a starters you don't have mean(est_sd) anywhere. The standard deviation for $y_i|\mu$ is given. Usually when people want to sample from the posterior they have one set of data $y$, but here you generate 1000 samples which causes some confusion. – ACE May 16 '16 at 8:36
• @Henry My bad, I forgot to include r. Fixed that. I had set mean = 5 as a way of setting the truth. I changed it to $\mu = 1$ to address Bjorn's concern – Greenparker May 16 '16 at 10:42
• @Björn Yes, $N(0,1)$ is the prior on $\mu$. I changed the true likelihood so that the prior is not awful. – Greenparker May 16 '16 at 10:45

I interpret your code as follows:

For $B = 1, 2, \dots 100$

(i) Generate $y \sim N(1_n, I_n)$

(ii) Generate $X \mid y \sim N(1_r [n+1]^{-1}n\bar{y}, [n+1]^{-1}I_r)$

(iii) Compute and store $\bar{X} = 1_r'X/r$

(i) and (ii) gives you independent samples from the joint distribution of $(y, X)$. Thus, \begin{align} \mathrm{Var}(\bar{X}) &= \mathbb E \mathrm{Var}(\bar{X} \mid y) + \mathrm{Var}(\mathbb E (\bar{X} \mid y))\\ & = \frac{1}{r(n+1)} + \mathrm{Var}\left( \bar{y}\frac{n}{n+1}\right) \\ &= \frac{1}{r(n+1)} + \frac{n}{(n+1)^2}. \end{align}

Thus, the variance of $\bar{X}$ is of the order $1 / (n + 1)$ so you would expect the sample standard deviation over the $B$ runs to estimate its square root fairly well, that is, to about $\sqrt{1/1001}$.

This also shows that the conjecture in the second question is false since, if the distributions of $\bar{X}$ and $X\mid y$ were the same, they would have the same variance.

• Thanks for the answer. I am unclear as to how you got $\frac{\bar{y}n}{n+1}$ as the expectation of $\bar{X}|y$. Shouldn't $E(\bar{X}|y)$ be $\bar{y}$? – Greenparker May 16 '16 at 20:52
• Your code says $X_i \sim N(\sum_i y_i / (n+1), \cdot) = N(\bar{y}n/(n+1), \cdot)$, and then $\mathbb E \bar{X} = \mathbb E X_i$. Ah, my last edit didn't save. Wait. Now it's probably coherent @Greenparker. – ekvall May 16 '16 at 20:56
• Oops, sorry, that was stupid of me. Ok as it stands, it makes sense to me. However, I feel like this is an oversimplified version of my problem, since I actually have $y|X = x \sim N(x1_n, I_n)$, and then the same $X|y$ as you. So my (i) and (ii) does not give the joint distribution of $(y,X)$ – Greenparker May 16 '16 at 21:16
• @Greenparker In that case, your code doesn't do what you want it to. Your code does what I've written here (I mean, barring any mistakes on my part). You are generating $y$ with fixed mean in every iteration. – ekvall May 16 '16 at 21:20
• You are right, the code isn't doing what I thought it was, and so your answer seems correct. Its possible that I have an inadequate understanding of the other question's data generation. I guess for the other question, the statement is false, because the variance for $\bar{X}$ is clearly not the same as the variance for $X$, so the distributions cannot be the same. Thanks! – Greenparker May 16 '16 at 21:24