I am trying to setup VB to do a weighted linear regression for vector observations. My setup is that I have $N$ numbers of $d$-dimensional vector observations. I would like to model the noise as being independent for an output observation but I would like each point to have its own variance. In view of that, my model is as follows:
$$ y_i \sim \mathcal{N}(T(x_i; \beta), \frac{\sigma^2}{w_i} {\textbf{I}}) \\ \beta \sim \mathcal{N} (\beta_{0}, \Sigma_{0}) \\ w_i \sim \mathcal{G} (a, b) \\ $$
Here $y_i$ is the one $d$-dimensional output observations and $x_i$ is the corresponding independent variable. They are related by some transformation $T$ which is parameterized by $\beta$. The output variance for a point $i$ is scaled by $w_i$ and $\sigma^2$ is some given global variance.
The full joint model can be written as:
$$ p(\beta, w, y |x) = \prod_{i=1}^{N} \bigg[p(y_i|x_i, w_i, H) \ p(w_i)\bigg] \ p(\beta) $$
Applying the mean field approximation $Q(w, \beta) = Q(w)Q(\beta)$ and getting the optimal distribution for $w_i$'s we need to take the expectation with respect to the $\beta$ parameters. We can work with the logarithm to make things a bit simpler
$$ \ln q(w^*) = \bigg\langle \ln \prod_{i=1}^{N} \bigg[p(y_i|x_i, w_i, \beta) \ p(w_i)\bigg] \ p(\beta)\bigg\rangle_{\beta} \\ =\bigg\langle \ln \prod_{i=1}^{N} \bigg[p(y_i|x_i, w_i, H) \ p(w_i)\bigg]\bigg\rangle_{\beta} + C $$
Now, we look at the term $p(y_i|x_i, w_i, H) \ p(w_i)$ in more detail. Expanding we have:
$$ \prod_{i=1}^{N} \frac{1}{\sqrt{2 \pi}} \big(\frac{\sigma^2}{w_i} {\textbf{I}}\big)^{-\frac{1}{2}} \exp \bigg(-0.5 \times [y_i - T(x_i)]^T \big(\frac{\sigma^2}{w_i} {\textbf{I}}\big)^{-1} \ [y_i - T(x_i)]\bigg) \ \bigg(\frac{b^a}{\Gamma(a)} w_i^{a-1} \exp{(-b w_i)}\bigg) $$
Dropping the constants we get:
$$ \propto \prod_{i=1}^{N} \big(\frac{\sigma^2}{w_i} {\textbf{I}}\big)^{-\frac{1}{2}} \exp \bigg(-0.5 \times [y_i - T(x_i)]^T \big(\frac{\sigma^2}{w_i} {\textbf{I}}\big)^{-1} \ [y_i - T(x_i)]\bigg) \bigg(w_i^{a-1} \exp{(-b w_i)}\bigg) \\ $$
Gathering some like terms we get:
$$ \propto \prod_{i=1}^{N} \big(\frac{\sigma^2}{w_i} {\textbf{I}}\big)^{-\frac{1}{2}} \ w_i^{a-1} \exp \bigg(-0.5 \times [y_i - T(x_i)]^T \big(\frac{\sigma^2}{w_i} {\textbf{I}}\big)^{-1} \ [y_i - T(x_i)] - b w_i\bigg) $$
Not sure if I should simplify this further and try and absorb the scalar $w_i$s into the matrix terms but it still does not help me somehow.
Now I have two issues. I was hoping at this point, this will start to look like some familiar distribution. It sort of looks like some Gamma distribution perhaps but I am not sure how I can simplify it to get it in that form.
My second issue is that after I have done that, I need to still take the expectation with respect to $Q(\beta)$ and I am not sure how to do that.
Any pointers on how to proceed next would be very appreciated!
[EDIT] So, I have been thinking about this. I am assuming I can write the joint model for each $i$ term as:
$$ p(y_i, w_i, \beta|x_i) = \prod_{i=1}^{N}p(y_i|x_i, w_i, \beta) p(w_i) \frac{p(\beta)}{N} $$
Now, for example, when computing the distribution for $w_i$ with the mean field approximation, I need to take the expectation wrt to $\beta$. So, expanding this, I have:
$$ \propto \sqrt{w_i} \exp{\bigg(-[y_i - T(x_i, \beta)]^T w_i [y_i - T(x_i, \beta)]\bigg)} \\ \bigg( w_i^{a-1} \exp{(-b w_i)}\bigg) \exp{\bigg(-(\beta - \beta_0)^T \Sigma_{ 0}(\beta - \beta_0) \bigg)} $$
Now, the first two exponentials come together nicely to make a Gamma distribution and I am not sure though how to manipulate the last exponential so that the posterior over $w_i$ gets some nice form.