# Metropolis-within-Gibbs for parametric inference of a regressive model

I have a regressive model of this form $$\begin{equation} y=f(\theta)+\varepsilon \end{equation}$$ to describe observations $$y$$, with noise $$\varepsilon$$ and a parametric function $$f$$ with parameters $$\theta$$. In reality the regressive model is much more complicated than what I've shown but that's not relevant for my question.

Now I am interested in the scenario where the parameters $$\theta$$ vary throughout the experiments, i.e. they have an inherent variability, and I assume that their distribution is multivariate normal $$\theta\sim\mathcal{N}(\mu,\Sigma)$$. Hence, I decided to build a hierarchical model in order to identify the hyperparameters $$\mu,\Sigma$$, having chosen a multivariate normal $$\mathcal{N}(\mu_0,\Sigma_0)$$ and an inverse Wishart $$\mathcal{W}^{-1}(\Psi_0,\nu_0)$$, respectively, as their hyperpriors.

Then, I have implemented a Metropolis-within-Gibbs routine in MATLAB, to sample these conditional posteriors: $$\begin{equation} \label{eq:20} \small \mu|\theta,\Sigma,y \sim \mathcal{N}\left((\Sigma_0^{-1} + n\Sigma^{-1})^{-1}( \Sigma_0^{-1}\mu_0 + n \Sigma^{-1} \theta ),(\Sigma_0^{-1} + n \Sigma^{-1})^{-1}\right) \end{equation}$$ for the mean vector, and $$\begin{equation} \label{eq:21} \small \Sigma|\theta,\mu, y\sim\mathcal{W}^{-1}(\Psi_0 +\sum_{i=1}^n (\theta_i - \mu) (\theta_i - \mu)^T,n+\nu_0) \end{equation}$$ for the covariance matrix, where $$n=1$$. The above choices are motivated by conjugate priors. During the iterative process of the Gibbs sampler, the Metropolis-Hastings is used to take one sample from the posterior distribution $$\theta|\mu,\Sigma,y$$ (which is also multivariate normal).

My current results show a very wide posterior for $$\Sigma$$ and I cannot understand why that is. I am still unsure of:

• how many degrees of freedom $$\nu_0$$ I should consider for the inverse Wishart;
• if it is correct to assume that $$n=1$$ since I am taking only one sample from the MH sampler;
• and if during the Gibbs routine I need to multiply the inverse Wishart sample by a factor $$(\nu_0-k)$$ as shown here ($$k$$ is the number of parameters $$\theta$$)?

I took 10000 samples to generate the below plots (one of them is a zoom of the variance of the parameter).

UPDATE According to this question my formulation might be fundamentally wrong, in the sense that I should be using a Normal Inverse Wishart instead of a separate MVN and inverse Wishart. I would like to confirm this, and gain further insight into how the Gibbs sampler needs to consider the conditional posteriors. Any answer to this specific question is very welcome.