Steps to figure out a posterior distribution when it might be simple enough to have an analytic form?

Question

This was also asked at Computational Science.

I am trying to compute a Bayesian estimate of some coefficients for an autoregression, with 11 data samples: $$ Y_{i} = \mu + \alpha\cdot{}Y_{i-1} + \epsilon_{i} $$ where $\epsilon_{i}$ is Gaussian with mean 0 and variance $\sigma_{e}^{2}$ The prior distribution on the vector $(\mu, \alpha)^{t}$ is Gaussian with mean $(0,0)$ and a diagonal covariance matrix with diagonal entries equal to $\sigma_{p}^{2}$.

Based on the autoregression formula, this means that the distribution of the data points (the $Y_{i}$) are normal with mean $\mu + \alpha\cdot{}Y_{i-1}$ and variance $\sigma_{e}^{2}$. Thus, the density for all of the data points $(Y)$ jointly (assuming independence, which is fine for the program I am writing), would be: $$ p(Y \quad | (\mu, \alpha)^{t}) = \prod_{i=2}^{11}\frac{1}{\sqrt{2\pi\sigma_{e}^{2}}}\exp{\frac{-(Y_{i} - \mu - \alpha\cdot{}Y_{i-1})^{2}}{2\sigma_{e}^{2}}}.$$

By Bayes' theorem, we can take the product of the above density with the prior density, and then we'll just need the normalizing constant. My hunch is that this should work out to be a Gaussian distribution, so we can worry about the normalizing constant at the end rather than explicitly calculating it with integrals over $\mu$ and $\alpha$.

This is the part I am having trouble with. How do I compute the multiplication of the prior density (which is multivariate) and this product of univariate data densities? The posterior needs to be purely a density of $\mu$ and $\alpha$, but I cannot see how you'll get that out of such a product.

Any pointers are really helpful, even if you just point me in the right direction and then I need to go and do the messy algebra (which is what I've already attempted several times).

As a starting point, here is the form of the numerator from Bayes' rule: $$ \frac{1}{(2\pi\sigma_{e}^{2})^{5}\cdot{}2\pi\sigma_{p}^{2}} \exp{\biggl [ \frac{1}{2\sigma_{e}^{2}}\sum_{i=2}^{11}(Y_{i} - \mu - \alpha\cdot{}Y_{i-1})^{2} - \frac{\mu^{2}}{2\sigma_{p}^{2}} - \frac{\alpha^{2}}{2\sigma_{p}^{2}} \biggr ] }.$$

The issue is how to see that this reduces down to a Gaussian density of $(\mu, \alpha)^{t}$.

Added

Ultimately, this boils down to the following general problem. If you are given some quadratic expression such as $$A\mu^{2} + B\mu\alpha + C\alpha^{2} + J\mu + K\alpha + L$$ how do you put that into a quadratic form $(\mu-\hat{\mu},\alpha-\hat{\alpha})Q(\mu-\hat{\mu},\alpha-\hat{\alpha})^{t}$ for some 2x2 matrix $Q$? It's simple enough in easy cases, but what process do you use to obtain the mean estimates, $\hat{\mu}$ and $\hat{\alpha}$?

Note, I tried the straightforward option of expanding the matrix formula and then trying to equate coefficients as above. The problem, in my case, is that the constant $L$ is zero, and then I end up getting three equations in two unknowns, so it is underdetermined to just match coefficients (even if I assume a symmetric quadratic form matrix).

Answer 1

The clue that was in my answer to the previous answer is to look at how I integrated out the parameters - because you will do exactly the same integrals here. You question assumes the variance parameters known, so they are constants. You only need to look at the $\alpha,\mu$ dependence on the numerator. To see this, note that we can write:

$$p(\mu,\alpha|Y)=\frac{p(\mu,\alpha)p(Y|\mu,\alpha)}{\int\int p(\mu,\alpha)p(Y|\mu,\alpha)d\mu d\alpha}$$ $$=\frac{\frac{1}{(2\pi\sigma_{e}^{2})^{5}\cdot{}2\pi\sigma_{p}^{2}} \exp{\biggl [ -\frac{1}{2\sigma_{e}^{2}}\sum_{i=2}^{11}(Y_{i} - \mu - \alpha\cdot{}Y_{i-1})^{2} - \frac{\mu^{2}}{2\sigma_{p}^{2}} - \frac{\alpha^{2}}{2\sigma_{p}^{2}} \biggr ] }}{\int\int \frac{1}{(2\pi\sigma_{e}^{2})^{5}\cdot{}2\pi\sigma_{p}^{2}} \exp{\biggl [ -\frac{1}{2\sigma_{e}^{2}}\sum_{i=2}^{11}(Y_{i} - \mu - \alpha\cdot{}Y_{i-1})^{2} - \frac{\mu^{2}}{2\sigma_{p}^{2}} - \frac{\alpha^{2}}{2\sigma_{p}^{2}} \biggr ] }d\mu d\alpha}$$

Notice how we can pull the first factor $\frac{1}{(2\pi\sigma_{e}^{2})^{5}\cdot{}2\pi\sigma_{p}^{2}}$ out of the double integral on the denominator, and it cancels with the numerator. We can also pull out the sum of squares $\exp{\biggl [ -\frac{1}{2\sigma_{e}^{2}}\sum_{i=2}^{11}Y_{i}^{2} \biggr ]}$ and it will also cancel. The integral we are left with is now (after expanding the squared term):

$$=\frac{\exp{\biggl [ -\frac{10\mu^2+\alpha^2\sum_{i=1}^{10}Y_{i}^{2}-2\mu\sum_{i=2}^{11}Y_i-2\alpha\sum_{i=2}^{11}Y_{i}Y_{i-1}+2\mu\alpha\sum_{i=1}^{10}Y_i}{2\sigma_{e}^{2}} - \frac{\mu^{2}}{2\sigma_{p}^{2}} - \frac{\alpha^{2}}{2\sigma_{p}^{2}} \biggr ] }}{\int\int \exp{\biggl [ -\frac{10\mu^2+\alpha^2\sum_{i=1}^{10}Y_{i}^{2}-2\mu\sum_{i=2}^{11}Y_i-2\alpha\sum_{i=2}^{11}Y_{i}Y_{i-1}+2\mu\alpha\sum_{i=1}^{10}Y_i}{2\sigma_{e}^{2}} - \frac{\mu^{2}}{2\sigma_{p}^{2}} - \frac{\alpha^{2}}{2\sigma_{p}^{2}} \biggr ] }d\mu d\alpha}$$

Now we can use a general result from the normal pdf.

$$\int \exp\left(-az^2+bz-c\right)dz=\sqrt{\frac{\pi}{a}}\exp\left(\frac{b^2}{4a}-c\right)$$ This follows from completing the square on $-az^2+bz$ and noting that $c$ does not depend on $z$. Note that the inner integral over $\mu$ is of this form with $a=\frac{10}{2\sigma^2_e}+\frac{1}{2\sigma^2_p}$ and $b=\frac{\sum_{i=2}^{11}Y_i-\alpha\sum_{i=1}^{10}Y_i}{\sigma_{e}^{2}}$ and $c=\frac{\alpha^2\sum_{i=1}^{10}Y_{i}^{2}-2\alpha\sum_{i=2}^{11}Y_{i}Y_{i-1}}{2\sigma_{e}^{2}}+ \frac{\alpha^{2}}{2\sigma_{p}^{2}}$. After doing this integral, you will find that the remaining integral over $\alpha$ is also of this form, so you can use this formula again, with a different $a,b,c$. Then you should be able to write your posterior in the form $\frac{1}{2\pi|V|^{\frac{1}{2}}}\exp\left[-\frac{1}{2}(\mu-\hat{\mu},\alpha-\hat{\alpha})V^{-1}(\mu-\hat{\mu},\alpha-\hat{\alpha})^T\right]$ where $V$ is a $2\times 2$ matrix

Let me know if you need more clues.

update

(note: correct formula, should be $10\mu^2$ instead of $\mu^2$)

if we look at the quadratic form you've written in the update, we notice there is $5$ coefficients ($L$ is irrelevant for posterior as we can always add any constant which will cancel in the denominator). We also have $5$ unknowns $\hat{\mu},\hat{\alpha},Q_{11},Q_{12}=Q_{21},Q_{22}$. Hence this is a "well posed" problem so long as the equations are linearly independent. If we expand the quadratic $(\mu-\hat{\mu},\alpha-\hat{\alpha})Q(\mu-\hat{\mu},\alpha-\hat{\alpha})^{t}$ we get:

$$Q_{11}(\mu-\hat{\mu})^2+Q_{22}(\alpha-\hat{\alpha})^2+2Q_{12}(\mu-\hat{\mu})(\alpha-\hat{\alpha})$$ $$=Q_{11}\mu^{2} + 2Q_{21}\mu\alpha + Q_{22}\alpha^{2} - (2Q_{11}\hat{\mu}+2Q_{12}\hat{\alpha})\mu - (2Q_{22}\hat{\alpha}+2Q_{12}\hat{\mu})\alpha +$$ $$+Q_{11}\hat{\mu}^2+Q_{22}\hat{\alpha}^2+2Q_{12}\hat{\mu}\hat{\alpha}$$

Comparing second order coefficient we get $A=Q_{11},B=2Q_{12},C=Q_{22}$ which tells us what the (inverse) covariance matrix looks like. Also we have two slightly more complicated equations for $\hat{\alpha},\hat{\mu}$ after substituting for $Q$. These can be written in matrix form as:

$$ -\begin{pmatrix}2A & B \\ B & 2C\end{pmatrix} \begin{pmatrix}\hat{\mu} \\ \hat{\alpha}\end{pmatrix} = \begin{pmatrix}J \\ K\end{pmatrix} $$

Thus the estimates are given by:

$$ \begin{pmatrix}\hat{\mu} \\ \hat{\alpha}\end{pmatrix} = -\begin{pmatrix}2A & B \\ B & 2C\end{pmatrix}^{-1}\begin{pmatrix}J \\ K\end{pmatrix}=\frac{1}{4AC-B^2}\begin{pmatrix}BK-2JC \\ BJ-2KA\end{pmatrix} $$

Showing that we do not have unique estimates unless $4AC\neq B^2$. Now we have: $$\begin{array}{c c} A=\frac{10}{2\sigma^2_e}+\frac{1}{2\sigma^2_p} & B=\frac{\sum_{i=1}^{10}Y_i}{\sigma_{e}^{2}} & C=\frac{\sum_{i=1}^{10}Y_{i}^{2}}{2\sigma^2_e}+\frac{1}{2\sigma^2_p} \\ J=-\frac{\sum_{i=2}^{11}Y_i}{\sigma_{e}^{2}} & K=-\frac{\sum_{i=2}^{11}Y_{i}Y_{i-1}}{\sigma_{e}^{2}} \end{array}$$

Note that if we define $X_i=Y_{i-1}$ for $i=2,\dots,11$ and take the limit $\sigma^2_p\to\infty$ then the estimates for $\mu,\alpha$ are given by the usual least squares estimate $\hat{\alpha}=\frac{\sum_{i=2}^{11}(Y_{i}-\overline{Y})(X_{i}-\overline{X})}{\sum_{i=2}^{11}(X_{i}-\overline{X})^2}$ and $\hat{\mu}=\overline{Y}-\hat{\alpha}\overline{X}$ where $\overline{Y}=\frac{1}{10}\sum_{i=2}^{11}Y_i$ and $\overline{X}=\frac{1}{10}\sum_{i=2}^{11}X_i=\frac{1}{10}\sum_{i=1}^{10}Y_i$. So the posterior estimates are a weighted average between the OLS estimates and the prior estimate $(0,0)$.