I need guidance on computing the posterior distribution of a hyperparameter in a specific multivariate normal model. Here's a brief description of my problem:
I have a dataset where the observed variable $\mathbf{A}_{N^2\times 1}$ follows a multivariate normal distribution $\mathcal{N}(\mathbf{A}|\boldsymbol{\mu}\;\mathbf{vec}(\mathbf{b}),\sigma_{\mathbf{A}}^2\mathbb{I})$, where $\boldsymbol{\mu}_{N^2\times K^2}$ is a latent variable and considered as another latent variable of the model. In this model, $\mathbf{vec}(\mathbf{b})$ is a latent variable which is basically a ${K^2\times 1}$ vector. The prior over $\mathbf{b}$ is a Gaussian with zero mean and covariance $\sigma_{\mathbf{b}}^2\mathbb{I}$. My goal is to compute the posterior distribution of the hyperparameter $\sigma_{\mathbf{b}}^2$.
In Bayesian statistics, the posterior distribution is proportional to the likelihood multiplied by the prior. However, I am facing challenges in integrating out the latent variable $\mathbf{b}$ to formulate the likelihood function $ p(\mathbf{A} | \sigma_{\mathbf{b}}^2)$ and in combining this with a prior distribution for $\sigma_{\mathbf{b}}^2$.
By integrating out $\mathbf{b}$ and using completing the square method:e $\mathcal{N}\Big(\mathbf{vec}(\mathbf{b})|\mu_{\mathbf{b}}, \Sigma_{\mathbf{b}}\Big)$ where $\Sigma_{\mathbf{b}}=\Big(\frac{1}{\sigma_{\mathbf{b}}^2}\mathbb{I}+\frac{1}{\sigma_{\mathbf{A}}^2}\boldsymbol{\mu}^T\boldsymbol{\mu}\Big)^{-1}$ and $\mu_{\mathbf{b}}=\frac{1}{\sigma_{\mathbf{A}}^2}\Sigma_{\mathbf{b}}^{-1}\boldsymbol{\mu}^T\mathbf{A}$.
My question is:
I am doing Gibbs Sampling, are there any recommended that could simplify the computation of this posterior distribution and would give me a closed-form distribution for posterior of $\sigma_{\mathbf{b}}^2$?
I derived that the posterior distribution should be:
$p(\sigma^2_{\mathbf{b}}|\mathbf{A}, \ldots) = \mathcal{N}\left(\mathbf{A} \middle| \boldsymbol{0}, \sigma_{\mathbf{A}}^2\mathbb{I}_{N^2} + \sigma_{\mathbf{b}}^2\boldsymbol{\mu}\boldsymbol{\mu}^T\right) \times \frac{d^c}{\Gamma(c)(\sigma^2_{\mathbf{b}})^{c+1}}\exp\left(-\frac{d}{\sigma^2_{\mathbf{b}}}\right).$
However, I find it non-intuitive how to proceed further since $\sigma^2_{\mathbf{b}}$ appears in the covariance matrix of a multivariate Gaussian distribution. I am seeking help on how to extract it in order to construct the posterior distribution of $\sigma^2_{\mathbf{b}}$
I would greatly appreciate any solutions, or examples of similar computations.
You already calculated your posterior distribution. More explicitly, the $\sigma_b^2$ dependence is captured by (the normalising prefactor is $\sigma_b^2$ independent): $$ p(\sigma_b^2|A,\mu,\sigma_A^2) \propto e^{-\frac{A^T\mu\mu^TA}{2\sigma_A^2\sigma_b^2}}\sqrt{|1+\sigma_b^2\mu\mu^T/\sigma_A^2|^{-1}}p_p(\sigma_b^2) $$ with $p_p(\sigma_b^2)$ the prior distribution, which I guess you chose to be an inverse gamma. In your case, this may not be judicious choice as it is not a conjugate prior to the posterior. It only works when $\sigma_A^2=0$ due to the determinant. You can make the determinant slightly more explicit: $$ \sqrt{|1+\sigma_b^2\mu\mu^T/\sigma_A^2|} = \prod_{i=1}^{N^2}\sqrt{|1+\sigma_b^2\lambda_i|} $$ with $\lambda_i$ the eigenvalues of $\mu\mu^T/\sigma_A^2$ (positive diagonalisable by the spectral theorem). Thus the conjugate prior should not have powers of $\sigma_b$ like for the gamma distribution but rather factors of this form.
Sampling $\sigma_b^2$ from this distribution would be tricky as it is less standard. Depending on what you originally want to sample, it would be best to just simulate normal distributions which define most of your variables.
Hope this helps.