\sim N \bigg( (\frac{\eta^{*T} \eta^*}{\sigma^2} + \frac{p(p + 1)}{2})^{-1} \eta^* \big(y - \eta \omega)\big)/\sigma^2,(\frac{\eta^{*T} \eta^*}{\sigma^2} + I_\frac{p(p + 1)}{2}/100)^{-1} \bigg)\end{aligned}$$\
where $\eta^*$ is a matrix containing the pairwise interactions of among the columns of $\eta$. Then set $\Omega = \frac{\Omega + \Omega^T}{2}$ *Step 4* Sample $\sigma^{-2}$ from a Gamma distribution: $$\begin{aligned}
\pi (\sigma^{-2} | \hbox{ --- } ) \sim Gamma \bigg( \frac{1 + n}{2}, \frac{1}{2} + \frac{1}{2} (y - \eta \omega - diag(\eta \Omega \eta^T))^T(y - \eta \omega - diag(\eta \Omega \eta^T)  \bigg)\end{aligned}$$\
*Step 5* Denote $\lambda_j$ the rows of $\Lambda$, for $j = 1,\cdots,p$. Sample $p$ conditionally independent posteriors: $$\begin{aligned}
\pi (\lambda_j | \hbox{ --- } ) \sim N \bigg( (D_j^{-1} + \frac{\eta^T\eta}{\sigma_j^2})^{-1} \eta^T \sigma_j^{-2} X^{(j)},(D_j^{-1} + \frac{\eta^T\eta}{\sigma_j^2})^{-1} \bigg)\end{aligned}$$\
where $X^{(j)} $ is the $j^{th}$ column of the matrix $X$, $D_j = diag(\tau_j^2 \psi_{j1} \phi_{j 1},\cdots, \tau_j^2 \psi_{jk} \phi_{j k})$. *Step 6* Sample $\psi_{jh}$ for $j = 1, \cdots, p$ and $h = 1 \cdots, k$ from independent Inverse Gaussian distributions: $$\begin{aligned}
\pi (\psi_{jh} | \hbox{ --- } ) \sim InvGauss \bigg( \tau_j \phi_{jh} , 1 \bigg)\end{aligned}$$\
*Step 7* Sample $\tau_{j}$ for $j = 1, \cdots, p$ from independent Generalized Inverse Gaussian distributions: $$\begin{aligned}
\pi (\tau_{j} | \hbox{ --- 