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])
#[1] 0.03159523
#[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$?

  • $\begingroup$ Your R code fails because object 'r' not found. It is also unclear why you have mean = 5 $\endgroup$
    – Henry
    May 16, 2016 at 6:37
  • $\begingroup$ 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... $\endgroup$
    – Björn
    May 16, 2016 at 7:52
  • $\begingroup$ 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. $\endgroup$
    – AnthonyC
    May 16, 2016 at 8:36
  • $\begingroup$ @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 $\endgroup$ May 16, 2016 at 10:42
  • $\begingroup$ @Björn Yes, $N(0,1)$ is the prior on $\mu$. I changed the true likelihood so that the prior is not awful. $\endgroup$ May 16, 2016 at 10:45

1 Answer 1


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.

  • $\begingroup$ 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}$? $\endgroup$ May 16, 2016 at 20:52
  • $\begingroup$ 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. $\endgroup$
    – KOE
    May 16, 2016 at 20:56
  • $\begingroup$ 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)$ $\endgroup$ May 16, 2016 at 21:16
  • 1
    $\begingroup$ @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. $\endgroup$
    – KOE
    May 16, 2016 at 21:20
  • 1
    $\begingroup$ 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! $\endgroup$ May 16, 2016 at 21:24

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.