Help with Bayesian derivation of normal model with conjugate prior

Question

I'm reading Marin's and Robert's Bayesian Essentials with R to work on the foundations of my understanding of Bayesian analysis. In the first chapter, about normal models, I'm having some trouble working the posterior distribution by myself.

For the model with known variance and conjugate prior for the mean $N(0, \sigma^2)$, the derivation seems straightforward (pages 29 and 30):

$$\pi(\mu\mid D_n) \propto \mathrm{exp}\{-\mu^2/2\sigma^2\}\mathrm{exp}\{-[n(\bar{x}-\mu)^2 + s^2]/2\sigma^2\}$$ $$\propto \mathrm{exp}\left\{-\frac{\mu^2 + n\mu^2 - 2n\mu\bar{x} +n\bar{x}^2 + s^2}{2\sigma^2}\right\}$$

$s^2$ drops off because it does not depend on the $\mu$ parameter; $n\bar{x}^2$ is almost what we need to complete the square, but not quite so, so I guess we also drop it off and substitute it for the proper term:

$$\propto \mathrm{exp}\left\{-\frac{(n+1)\mu^2 - 2n\mu\bar{x}}{2\sigma^2}\right\}$$ $$\propto \mathrm{exp}\left\{-\frac{(n+1)\left(\mu^2 - \frac{2n\mu\bar{x}}{n+1} + \frac{n^2\bar{x}^2}{(n+1)^2}\right)}{2\sigma^2}\right\}$$ $$\propto \mathrm{exp}\left\{-\frac{(n+1)\left(\mu - \frac{n\bar{x}}{n+1}\right)^2}{2\sigma^2}\right\}$$

So far, so good. But when it comes to the unknown mean, unknown variance model with prior $N(0, \sigma^2)$ and $\frac{1}{\sigma^2} \sim Exp(1)$, we can't drop off the $n\bar{x}^2$ term anymore because now the denominator inside the exponential function depends on the parameter $\sigma^2$.

So I have come to a similar point deriving the joint posterior for $\mu$ and $\sigma^2$:

$$\propto \mathrm{exp}\left\{-\frac{2 + \mu^2 + n\mu^2 - 2n\mu\bar{x} +n\bar{x}^2 + s^2}{2\sigma^2}\right\}\sigma^{-(n+5)}$$ $$\propto \sigma^{-1/2}\mathrm{exp}\left\{-\frac{(n+1)\mu^2 - 2n\mu\bar{x} +n\bar{x}^2}{2\sigma^2}\right\}\mathrm{exp}\left\{-\frac{2+s^2}{2\sigma^2}\right\}(\sigma^2)^{-(\frac{n+2}{2})-1}$$

The part to the right is already the core of the Inverse-Gamma function derived in the book, but I can't see how to complete the square for the normal density function without adding another term to the second Inverse-Gamma parameter. Indeed, in the general case, this parameter value should also have a $\frac{n\bar{x}^2}{n+1}$ term (as, e.g., https://en.wikipedia.org/wiki/Normal_distribution#With_unknown_mean_and_unknown_variance).

So I decided to compare both solutions by simulation using Metropolis-Hastings to obtain posterior draws of both parameters. Surprisingly (or not), the book solution uses the correct formula for the scale parameter for the Inverse-Gamma, $\frac{2 + s^2}{2}$, and if we add the term based on the sample mean, the density curve does not quite fit (Edit: not really, see below) . What am I missing here?

---

Edit: I'm starting to think that the book might indeed be mistaken. Completing the square in the first exponential makes it necessary to add a $\frac{n\bar{x}^2}{n+1}$ term that can only go in the second exponential term. So I tried my hand again with simulations, adding a JAGS model with the same definition to compare to my homemade solution (which had a coding error when I posted the question), and now it seems that the correct formula for the second Inverse-Gamma parameter is indeed $\frac{2 + s^2 + \frac{n\bar{x}^2}{n+1}}{2}$. Here's a plot with the simulated and theoretical densities for the posterior distribution (data = rnorm(100, 50, 10)):

Answer 1

[Here is the full text from our book]

The case of a normal distribution with a known variance being quite unrealistic, we now consider the general case of an iid sample $$\mathscr{D}_n=(x_1,\ldots,x_n)$$ from the normal distribution $\mathscr{N}(\mu,\sigma^2)$ and $\theta=(\mu,\sigma^2)$. Keeping the same prior distribution $$\mathscr{N}\left(0,\sigma^{2}\right)$$ on $\mu$, which then appears as a conditional distribution of $\mu$ given $\sigma^2$, {\em i.e.}, relies on the generic decomposition $$ \pi(\mu,\sigma^2) = \pi(\mu|\sigma^2) \pi(\sigma^2)\,, $$ we have to introduce a further prior distribution on $\sigma^2$. To make computations simple at this early stage, we choose an exponential $\mathscr{E}(1)$ distribution on $\sigma^{-2}$. This means that the random variable $\omega=\sigma^{-2}$ is distributed from an exponential $\mathscr{E}(1)$ distribution, the distribution on $\sigma^2$ being derived by the usual change of variable technique, $$ \pi(\sigma^2) = \exp(-\sigma^{-2})\,\left|\dfrac{\text{d}\sigma^{-2}}{\text{d}\sigma^2}\right| = \exp(-\sigma^{-2})\,(\sigma^2)^{-2}\,. $$ (This distribution is a special case of an inverse gamma distribution, namely $\mathcal{IG}(1,1)$.) The corresponding posterior density on $\theta$ is then given by \begin{eqnarray*}\label{eq:conjunor} \pi((\mu,\sigma^2)|\mathscr{D}_n) &\propto& \pi(\sigma^2)\times\pi(\mu|\sigma^2)\times f(\mathscr{D}_n| \mu,\sigma^2)\nonumber\\ & \propto & (\sigma^{-2})^{1/2+2}\, \exp\left\{-(\mu^2 + 2)/2\sigma^2\right\}\nonumber\\ && \times (\sigma^{-2})^{n/2}\,\exp \left\{-\left(n(\mu-\overline{x})^2 + s^2 \right)/2\sigma^2\right\} \\ &\propto& (\sigma^2)^{-(n+5)/2}\exp\left\{-\left[(n+1) (\mu-n\bar x/(n+1))^2+(2+s^2)\right]/2\sigma^2\right\}\nonumber\\ &\propto& (\sigma^2)^{-1/2}\exp\left\{-(n+1)[\mu-n\bar x/(n+1)]^2/2\sigma^2\right\}\,.\nonumber\\ &&\times (\sigma^2)^{-(n+2)/2-1}\exp\left\{-(2+s^2)/2\sigma^2\right\}\,.\nonumber \end{eqnarray*} Therefore, the posterior on $\theta$ can be decomposed as the product of an inverse gamma distribution on $\sigma^2$, $$\mathscr{IG}((n+2)/2,[2+s^2]/2)$$, which is the distribution of the inverse of a gamma $$\mathscr{G}((n+2)/2,[2+s^2]/2)$$ random variable, and, conditionally on $\sigma^2$, a normal distribution on $\mu$, $$\mathscr{N} (n\bar x/(n+1),\sigma^2/(n+1)).$$ The interpretation of this posterior is quite similar to the case when $\sigma$ is known, with the difference that the variability in $\sigma$ induces more variability in $\mu$, the marginal posterior in $\mu$ being then a Student's $t$ distribution $$ \mu|\mathscr{D}_n \sim \mathscr{T}\left(n+2,n\bar x/(n+1),(2+s^2)/(n+1)(n+2)\right)\,, $$ with $n+2$ degrees of freedom, a location parameter proportional to $\bar x$ and a scale parameter (almost) proportional to $s$.

But then, indeed, there is a mistake in the derivation of the marginal on $\sigma^2$ which, as you point out, should include a discrepancy between the prior mean ($0$) and the sample mean. The main derivation should thus be \begin{eqnarray*}\label{eq:conjucor} \pi((\mu,\sigma^2)|\mathscr{D}_n) &\propto& \pi(\sigma^2)\times\pi(\mu|\sigma^2)\times f(\mathscr{D}_n| \mu,\sigma^2)\nonumber\\ & \propto & (\sigma^{-2})^{1/2+2}\, \exp\left\{-(\mu^2 + 2)/2\sigma^2\right\}\nonumber\\ && \times (\sigma^{-2})^{n/2}\,\exp \left\{-\left(n(\mu-\overline{x})^2 + s^2 \right)/2\sigma^2\right\} \\ &\propto& (\sigma^2)^{-(n+5)/2}\exp\left\{-\left[(n+1) (\mu-n\bar x/(n+1))^2+(2+s^2)\right]/2\sigma^2\right\}\nonumber\\ &\propto& (\sigma^2)^{-1/2}\exp\left\{-(n+1)[\mu-n\bar x/(n+1)]^2/2\sigma^2\right\}\,.\nonumber\\ &&\times (\sigma^2)^{-(n+2)/2-1}\exp\left\{-\left(2+s^2+\frac{n}{n+1}{\bar x}^2\right)/2\sigma^2\right\}\,.\nonumber \end{eqnarray*} This means that the posterior on $\sigma^2$ is an inverse gamma distribution

$$\mathscr{IG}\left((n+2)/2,\left[2+s^2+\frac{n}{n+1}{\bar x}^2\right]/2\right)$$

Mea culpa! And apologies...