# Help with Bayesian derivation of normal model with conjugate prior

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)): • A good place to begin would be to get the algebra right: I'm afraid that no step of the derivation you lay out is correct. It might be a little easier to carry out in terms of logarithms so that you don't have to keep writing everything as an exponential.
– whuber
Sep 21, 2017 at 21:03
• whuber, do you mean that even for the first part (derivation of the posterior with known variance), the algebra is wrong? I noticed I made some mistakes in the Latex code that were repeated thanks to copy and paste, so maybe it makes more sense now. Sep 21, 2017 at 23:25

[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...

• Professor Robert, thank you very much for your input! I understand my problem is of algebraic nature: how do we complete the square to go from the second to the third line in the definition of the joint posterior in your book. Using the general case to work this specific case out ($\mu \sim N(0, \sigma^2)$ and $\sigma^2 \sim IG(1, 1)$), and using the notation from Wikipedia, we have: $n_0 = 1$, $\nu_0=2$, $\nu_0 \sigma_0 = 2$ and $\mu_0 = 0$. So, shouldn't we obtain an Inverse Gamma with parameters $\frac{\nu_0 + n}{2} = \frac{2 + n}{2}$ (same as the book) and ... Sep 24, 2017 at 16:21
• ... $\frac{\nu_0\sigma_0 +s^2 + \frac{n_0n}{n_0 +n} (\mu_0 - \bar{x})^2}{2} = \frac{2 +s^2 + \frac{n}{1 +n} \bar{x}^2}{2}$ -- containing the additional term necessary to complete the square but not present in the book derivation? (On a side note, I am big fan of your books and blog!) Sep 24, 2017 at 16:23
• You are completely right! Sep 24, 2017 at 17:11