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

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}$.

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

-
My answer to [this question]( stats.stackexchange.com/questions/22852/…) may be helpful. Note that you need a prior for your first observation - the iterations stop there. –  probabilityislogic Feb 18 '12 at 3:33
I don't see why I need it in this case. I'm supposed to treat the time intervals like they are conditionally independent given the observation. Notice that the product of the joint density is just from $i=2..11$. I don't think I'm supposed to get a sequentially updated formula here, just a single formula for the posterior $p((\mu,\alpha)^{t}\quad |Y)$. –  EMS Feb 18 '12 at 3:42
The "multivariate" in the prior $p(\alpha,\mu)$ is not in contradiction with the "univariate" in the data densities, because they are densities in the $y_i$'s. –  Xi'an Feb 18 '12 at 9:02
Thanks for updating your question. I will update my answer accordingly as i properly understand where you're getting stuck. –  probabilityislogic Feb 20 '12 at 22:22
Thanks. I appreciate your patience with me; I know it is elementary but your comments are quite helpful. I'm looking forward to any comments on the quadratic form part of the problem. –  EMS Feb 21 '12 at 1:40

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)$.

-
This isn't particularly helpful because I mentioned specifically that it's not the denominator that matters here. The denominator is just a normalizing constant, which will be obvious once you reduce the numerator to a Gaussian form. So tricks for evaluating the integrals in the denominator are mathematically really cool, but just not needed for my application. The only issue I need resolution with is manipulating the numerator. –  EMS Feb 18 '12 at 8:04
This answer gives you both numerator and denominator. The numerator exhibits the proper second degree polynomial in $(\alpha,\mu)$ that leads to the normal quadratic form, as stressed by probabilityislogic. –  Xi'an Feb 18 '12 at 9:05
@ems - by calculating the normalising constant you will construct the quadratic form required. it will contain the terms needed to compllete the square –  probabilityislogic Feb 18 '12 at 9:14
I don't understand how this gives you the quadratic form. I've worked out the two integrals in the denominator using the Gaussian integral identity that you posted. In the end, I just get a huge, messy constant. There doesn't seem to be any clear way to take that constant and turn it into something times a determinant to the 1/2 power, etc. Not to mention I don't see how any of this explains how to calculate the new 'mean vector' $(\hat{\mu},\hat{\alpha})^{t}$.. This is what I was asking for help for in the original question. –  EMS Feb 20 '12 at 19:39
Thanks tremendously for the detailed addition. I was making some silly errors when trying to do the algebra to figure out the quadratic form. Your comments about the relation to the OLS estimator are highly interesting and appreciated as well. I think this will speed up my code because I'll be able to draw samples from an analytic form that has built-in, optimized methods. My original plan was to use Metropolis-Hastings to sample from this, but it was very slow. Thanks! –  EMS Feb 21 '12 at 8:59