Multi-dimensional Integration Let $C=\sum_{s=1}^tx_sx_s'+aD^{-1}$, which is symmetric positive definite, vector $x,w\in\mathbb{R}^n$ and a scalar $y\in\mathbb{R}$. The integral I am trying to solve is as follows:
$$\sqrt{\frac{a\sqrt{\sigma^2}}{2\pi}}\frac{\int_{\mathbb{R}^n}e^{-\frac{a\sqrt{\sigma^2}}{2}(w' x_{t+1}-y)^2}e^{-\frac{a\sqrt{\sigma^2}}{2}\sum_{s}(w'x_s-y_s)^2-\frac{a^2\sqrt{{\sigma^2}}}{2} w' D^{-1}w}\mathrm{d}w}{\int_{\mathbb{R}^n}e^{-\frac{a\sqrt{\sigma^2}}{2}\sum_{s}(w'x_s-y_s)^2-\frac{a^2\sqrt{{\sigma^2}}}{2} w' D^{-1}w}\mathrm{d}w}$$
By using the arguments mentioned here I can solve the integral in the denominator by letting $b=2\sum_{s=1}^tx_sy_s$ and $\eta=\frac{a\sqrt{\sigma^2}}{2}$:
$\int_{\mathbb{R}^n}e^{-\eta w'Cw+\eta b'w-\eta\sum_{s}y_s^2}\mathrm{d}w=e^{-\eta\sum_{s}y_s^2}\prod_{i=1}^n\int_{\mathbb{R}}e^{-\eta\mu^i(w^{i})^2+\eta w^i(Pb)^i}\mathrm{d}w^i=e^{-\eta \sum_{s}y_s^2}\prod_{i=1}^n\left(\frac{\sqrt{\pi}e^{(\eta b^i)^2/4\eta\mu^i}}{\sqrt{\eta\mu^i}}\right)=e^{-\eta \sum_{s}y_s^2}\frac{\pi^{n/2}e^{\frac{\eta^2 b'b}{4\eta (\mu^1+...+\mu^n) }}}{\sqrt{\det \eta C}}$
To solve the numerator I write:
$$\int_{\mathbb{R}^n} e^{-\eta\left(\sum_sy_s-y^2\right)+\eta\left(\sum_s y_s x_s- yx_{t+1}\right)'w-\eta w'\left(aD^{-1}+\sum_s x_tx_t'+x_{t+1}x_{t+1}'\right)w}\mathrm{d}w
$$$$=\int_{\mathbb{R}^n} e^{-\eta\left(\sum_sy_s-y^2\right)+2\eta\left(\sum_s y_s x_s- yx_{t+1}\right)'w-\eta w'\left(aD^{-1}+\sum_s^{t+1} x_tx_t'\right)w}\mathrm{d}w$$
By using similar arguments as used to solve the denominator I write:
$$e^{-\eta \sum_{s}y_s^2-\eta y^2}\frac{\pi^{n/2}e^{\frac{\eta^2 b++'b++}{4\eta (\mu^1+...+\mu^n+x^1+...+x^n)}}}{\sqrt{\det 4\eta C++}}$$
So, I am left with the following:
$\sqrt{\frac{\eta}{\pi}}\left(e^{-\eta \sum_{s}y_s^2-\eta  y^2}\frac{\pi^{n/2}e^{\frac{\eta^2 b++'b++}{4\eta (\mu^1+...+\mu^n+x^1+...+x^n)}}}{\sqrt{\det C++}}\div e^{-\eta \sum_{s}y_s^2}\frac{\pi^{n/2}e^{\frac{\eta^2 b'b}{4\eta (\mu^1+...+\mu^n)}}}{\sqrt{\det 4\eta C}}\right)=\sqrt{\frac{\eta}{\pi}}\frac{e^{\frac{\eta y^2+(\eta^2 b++'b++)}{{4\eta (\mu^1+...+\mu^n+x^1+...+x^n)}}-\frac{\eta^2 b'b}{{4\eta (\mu^1+...+\mu^n)}}}\sqrt{\det 4\eta C}}{\sqrt{\det 4\eta C++}}$
where $C++=aD^{-1}+\sum_s^{t+1} x_tx_t'$ and $b++=2\sum_s y_s x_s- 2yx_{t+1}$. So I write:
$$\sqrt{\frac{\eta}{\pi}}\frac{e^{\frac{\eta y^2+(\eta^2 b++'b++)}{{4\eta (\mu^1+...+\mu^n+x^1+...+x^n)}}-\frac{\eta^2 b'b}{{4\eta (\mu^1+...+\mu^n)}}}\sqrt{\det 4\eta C}}{\sqrt{\det 4\eta C++}}=\sqrt{\frac{\eta}{\pi}}\left(\frac{e^{\frac{\eta y^2 +\eta^2(b-2yx_{t+1})'(b-2yx_{t+1})}{4\eta (\mu^1+...+\mu^n+x^1+...+x^n)}-\frac{\eta^2b'b}{4\eta (\mu^1+...+\mu^n)}}\sqrt{\det C}}{\sqrt{\det (C+x_{t+1}x_{t+1}')}}\right)$$

Maybe I can proceed further by using the information given here and (assuming $x_{t+1}x_{t+1}'$ is positive definite(Ideally I would like to avoid this)) write:
  $$\sqrt{\frac{\eta}{\pi}}\left(\frac{e^{\frac{\eta y^2 +\eta^2(b-2yx_{t+1})'(b-2yx_{t+1})}{4\eta (\mu^1+...+\mu^n+x^1+...+x^n)}-\frac{\eta^2b'b}{4\eta (\mu^1+...+\mu^n)}}\sqrt{\det C}}{\sqrt{\det (C+x_{t+1}x_{t+1}')}}\right)=\sqrt{\frac{\eta}{\pi}}\left(\frac{e^{\frac{\eta y^2 +\eta^2(b-2yx_{t+1})'(b-2yx_{t+1})}{4\eta (\mu^1+...+\mu^n+x^1+...+x^n)}-\frac{\eta^2b'b}{4\eta (\mu^1+...+\mu^n+x^1+...+x^n)}}\sqrt{\det C}}{\sqrt{(4\eta)^n\det C\det(I+L'C^{-1}L)}}\right)=\sqrt{\frac{\eta}{\pi}}\left(\frac{e^{\frac{\eta y^2 +\eta^2(b-2yx_{t+1})'(b-2yx_{t+1})}{4\eta (\mu^1+...+\mu^n+x^1+...+x^n)}-\frac{\eta^2b'b}{4\eta (\mu^1+...+\mu^n)}}}{\sqrt{(4\eta)^n\det(I+L'C^{-1}L)}}\right)=\sqrt{\frac{\eta}{\pi}}\left(\frac{e^{\frac{\eta y^2 +\eta^2(b-2yx_{t+1})'(b-2yx_{t+1})-\eta^2b'b4\eta (x^1+...+x^n)}{4\eta (\mu^1+...+\mu^n+x^1+...+x^n)}}}{\sqrt{(4\eta)^n\det(I+L'C^{-1}L)}}\right)$$

Can someone please tell me: what I have done so far is correct and how can I proceed to simplify whats left? 
Remark: I am hoping to get a Gaussian density. Once I solve, then I can analytically get the mean and variance. It may be worth noting that $C+x_{t+1}x_{t+1}$ is also symmetric and positive definite. 
 A: To obtain the predictive distribution for Normal/Gaussian likelihood with sequence $S$ we need to solve the following:
\begin{equation}\label{pred}
 p(y|x_T,S_{T-1})=\int_{\mathbb{R}^n}p(y|x_T,w)p(w|S_{T-1})\mathrm{d}w
\end{equation}
where the (uniformly initialised) prior distribution is:
\begin{equation}\label{prior}
    p(w)=\left(\frac{a}{4\sigma^2}\right)^n\exp\left(-\frac{a}{2\sigma^2}w'D^{-1}w\right)
\end{equation}
and the posterior is:
\begin{equation}\label{post}
  p(w|S_{T-1})=  \frac{\left(\prod_{t=1}^{T-1}p(y_t|x_t,w)\right)p(w)}{\int_{\mathbb{R}^n}\left(\prod_{t=1}^{T-1}p(y_t|x_t,w)\right)p(w)\mathrm{d}w}
\end{equation}
Thus, the predictive distribution at time $T$ for $y$ given the sequence $S_{T-1}=x_1,y_1,..,x_{T-1},y_{T-1}$ requires evaluation of the following integral:
\begin{equation}
\frac{\int_{\mathbb{R}^n}\frac{1}{\sqrt{2\pi\sigma^2}}e^{\frac{(w'x_T-y)^2}{2\sigma^2}}\prod_{t=1}^{T-1}\frac{1}{\sqrt{2\pi\sigma^2}}e^{\frac{(w'x_t-y_t)^2}{2\sigma^2}}\exp\left(-\frac{a}{2\sigma^2}w'D^{-1}w\right)\mathrm{d}w}{\int_{\mathbb{R}^n}\prod_{t=1}^{T-1}\frac{1}{\sqrt{2\pi\sigma^2}}e^{\frac{(w'x_t-y_t)^2}{2\sigma^2}}\exp\left(-\frac{a}{2\sigma^2}w'D^{-1}w\right)\mathrm{d}w}
\end{equation}
Let $\eta=\frac{1}{2\sigma^2}$ and,
\begin{equation}\label{id1}
L_T^w=  \sum_{t=1}^{T}(y_t-w'x_t)^2=\sum_{t=1}^{T}y_t^2-2w'\left(\sum_{t=1}^{T}x_ty_t\right)+w'\left(\sum_{t=1}^{T}x_tx_t'\right)w
\end{equation}
For any $x,b\in\mathbb{R}^n$ and a symmetric positive definite matrix $A$:
$$x'Ax-2b'x=(x-A^{-1}b)'A(x-A^{-1}b)-b'A^{-1}b$$
Expanding quadratic form:
$$
(x-A^{-1}b)'A(x-A^{-1}b)=x'Ax-2b'A^{-1}Ax+b'A^{-1}AA^{-1}b
$$
$$=x'Ax-2b'x+b'A^{-1}b$$
immediately follows:
\begin{multline}\label{id2}
    w'\left(\sum_{t=1}^{T-1}x_tx_t'+aD^{-1}\right)w-2w'\left(\sum_{t=1}^{T-1}x_ty_t\right)=\\\left(w-\left(\sum_{t=1}^{T-1}x_ty_t\right)'\left(\sum_{t=1}^{T-1}x_tx_t'+aD^{-1}\right)^{-1}\right)'\left(\sum_{t=1}^{T-1}x_tx_t'+aD^{-1}\right)\\\left(w-\left(\sum_{t=1}^{T-1}x_ty_t\right)'\left(\sum_{t=1}^{T-1}x_tx_t'+aD^{-1}\right)^{-1}\right)-\\\left(\sum_{t=1}^{T-1}x_ty_t\right)'\left(\sum_{t=1}^{T-1}x_tx_t'+aD^{-1}\right)^{-1}\left(\sum_{t=1}^{T-1}x_ty_t\right)
\end{multline}
If an algorithm follows Bayesian strategy with Gaussian likelihood and prior such that absolute value of the each element of the weight vector is not zero, $w_0$ is initialised uniformly and $a > 0$, then the posterior distribution is:
$$
\mathcal{N}\left(\left(\sum_{t=1}^{T-1}x_t y_t\right)'\left(\sum_{t=1}^{T-1}x_t x_t'+aD^{-1}\right)^{-1},\frac{1}{2\sigma^2}\left(\sum_{t=1}^{T-1}x_t x_t'+aD^{-1}\right)^{-1}\right)
$$
Expanding posterior ignoring the normalising constant we get:
$$
p(w|S_{T-1})\propto \exp \left(-\eta\sum_{t=1}^{T-1}(y_t-w'x_t)^2-a\eta w'D^{-1} w\right)
$$
\begin{multline*}
=\exp\left(-\eta\Bigg(w'\left(\sum_{t=1}^{T-1}x_tx_t'+aD^{-1}\right)w-2w'\sum_{t=1}^{T-1}x_ty_t+\sum_{t=1}^{T-1}y_t^2\right)\Bigg)
\end{multline*}
\begin{multline*}
 =   \exp\left(-\eta\Bigg(w-\left(\sum_{t=1}^{T-1}x_ty_t\right)'\left(\sum_{t=1}^{T-1}x_tx_t'+aD^{-1}\right)^{-1}\right)'\\\left(\sum_{t=1}^{T-1}x_tx_t'+aD^{-1}\right)\left(w-\left(\sum_{t=1}^{T-1}x_ty_t\right)'\left(\sum_{t=1}^{T-1}x_tx_t'+aD^{-1}\right)^{-1}\right)\\-\left(\sum_{t=1}^{T-1}x_ty_t\right)'\left(\sum_{t=1}^{T-1}x_tx_t'+aD^{-1}\right)^{-1}\left(\sum_{t=1}^{T-1}x_ty_t\right)\\+\left(\sum_{t=1}^{T-1}x_ty_t\right)'\left(\sum_{t=1}^{T-1}x_tx_t'+aD^{-1}\right)^{-1}\left(\sum_{t=1}^{T-1}x_ty_t\right)+\sum_{t=1}^{T-1}y_t^2\Bigg)
\end{multline*}
\begin{multline}
\propto\exp\left(-\eta\Bigg(w-\left(\sum_{t=1}^{T-1}x_ty_t\right)'\left(\sum_{t=1}^{T-1}x_tx_t'+aD^{-1}\right)^{-1}\right)'\\\left(\sum_{t=1}^{T-1}x_tx_t'+aD^{-1}\right)\left(w-\left(\sum_{t=1}^{T-1}x_ty_t\right)'\left(\sum_{t=1}^{T-1}x_tx_t'+aD^{-1}\right)^{-1}\right)
\end{multline}
The last proportionality can be recognised as probability density function of the multivariate Normal distribution.
This posterior distribution can be thought of online variant of the posterior obtained by . Since the posterior predictive distribution is a weighted average over parameter space where each parameter is weighted by its posterior probability, thus the predictive distribution is as follows:
$$
    \mathcal{N}\Bigg(\left(\sum_{t=1}^{T-1}x_t y_t\right)'\left(\sum_{t=1}^{T-1}x_t x_t'+aD^{-1}\right)^{-1}x_T,\frac{1}{2\sigma^2}x_T\left(\sum_{t=1}^{T-1}x_t x_t'+aD^{-1}\right)^{-1}x_T\Bigg)
$$
