Consider a multivariate vairable regression model, $$Y=XB+E,$$ where $Y=(y_1,y_2,\dots ,y_n{)}^{\top }$ is the $n\times q$ response matrix, $X=(x_1,\dots ,x_p)$ is the $n\times p$ design matrix, $B=({\beta }_1,\dots ,{\beta }_p{)}^{\top }$ is a $p\times q$ unknown coefficient matrix, and $E=(e_1,e_2,\dots ,e_n{)}^{\top }$ is an $n\times q$ error matrix with $e_i\stackrel{i.i.d}{\sim }{N}_q(0,\Sigma )$. To identify the variables selected among $p$ regressors, (Brown and Vannucci 1998) introduced a latent binary vector $\xi =({\xi }_1,\dots ,{\xi }_p{)}^{\top }$ with $$\begin{array}{c}{\xi }_j=\{\begin{array}{cc}1& \text{if}{\beta }_j\text{is active}\\ 0& \text{if}{\beta }_j\text{is not active}\end{array}.\end{array}$$ The conjugate priors are considered as follows: $$\begin{array}{ccc}& & \pi (B_{\xi }|\Sigma ,\xi )\sim MN(0,H_{\xi },\Sigma )\\ & & \pi \left(\Sigma \right)\sim IW(Q,\delta )\\ & & \pi \left(\xi \right)\stackrel{ind}{\sim }\prod _{j=1}^p\text{Ber}\left({\omega }_j\right),\end{array}$$ where $Q$, $H_{\xi }$ and $\delta$ are hyperparameters. Then, by Bayesian theorem, the marginal posterior distribution can be obtained by $$\begin{array}{ccc}& & m\left(\xi |Y\right)\propto m\left(Y|\xi \right)\pi \left(\xi \right)\\ & =& \iint f(Y|B_{\xi },\Sigma )\pi (B_{\xi }|\Sigma ,\xi )\pi \left(\Sigma \right)dBd\Sigma \pi \left(\xi \right)\\ & =& \iint (2\pi {)}^{-\frac{nq}{2}}|\Sigma {|}^{-\frac{n}{2}}\exp [-\frac{1}{2}\text{t}r\{(Y-X_{\xi }{B}_{\xi }{)}^{\top }(Y-X_{\xi }{B}_{\xi }\left)\right\}{\Sigma }^{-1}]\times \\ & & (2\pi {)}^{-\frac{|\xi |q}{2}}|\Sigma {|}^{-\frac{|\xi |}{2}}|H_{\xi }{|}^{-\frac{q}{2}}\exp \left[\frac{1}{2}\text{t}r\right\{B_{\xi }^{\top }{H}_{\xi }^{-1}{B}_{\xi }{\Sigma }^{-1}\left\}\right]\times \\ & & |Q{|}^{\frac{\delta }{2}}{2}^{-\frac{1}{2}\delta q}{\xi }_q^{-1}\left(\frac{\delta }{2}\right)|\Sigma {|}^{-\frac{\delta +q+1}{2}}\exp \{-\frac{1}{2}\text{t}r(Q{\Sigma }^{-1}\left)\right\}d{B}_{\xi }d\Sigma \pi \left(\xi \right)\\ & \propto & (2\pi {)}^{-\frac{|\xi |q}{2}}{|H_{\xi }|}^{-\frac{q}{2}}\iint |\Sigma {|}^{-\frac{1}{2}(n+|\xi |+\delta +q+1)}\exp \{-\frac{1}{2}\text{t}r\left(Q{\Sigma }^{-1}\right)\}\\ & & \exp [-\frac{1}{2}\text{t}r\left\{\left(B_{\xi }^{\top }\right(X_{\xi }^{\top }{X}_{\xi }+H_{\xi }^{-1})B_{\xi }-2B_{\xi }^{\top }{X}_{\xi }^{\top }Y){\Sigma }^{-1}\right\}]\times \\ & & \exp [-\frac{1}{2}\text{t}r\left\{Y^{\top }Y{\Sigma }^{-1}\right\}]d{B}_{\xi }d\Sigma \pi \left(\xi \right)\\ & =& (2\pi {)}^{-\frac{|\xi |q}{2}}{|H_{\xi }|}^{-\frac{q}{2}}\iint |\Sigma {|}^{-\frac{1}{2}(n+|\xi |+\delta +q+1)}\exp \{-\frac{1}{2}\text{t}r\left[\right(Y^{\top }Y+Q\left){\Sigma }^{-1}\right]\}\times \\ & & \exp [-\frac{1}{2}\text{t}r\left\{(B_{\xi }^{\top }{K}_{\xi }{B}_{\xi }-2B_{\xi }^{\top }{K}_{\xi }{K}_{\xi }^{-1}{X}_{\xi }^{\top }Y){\Sigma }^{-1}\right\}]d{B}_{\xi }d\Sigma \pi \left(\xi \right)\\ & =& (2\pi {)}^{-\frac{|\xi |q}{2}}{|H_{\xi }|}^{-\frac{q}{2}}\iint |\Sigma {|}^{-\frac{1}{2}(n+|\xi |+\delta +q+1)}\times \\ & & \exp \{-\frac{1}{2}\text{t}r[(Y^{\top }Y-{\tilde{B}}_{\xi }^{\top }{K}_{\xi }{\tilde{B}}_{\xi }+Q){\Sigma }^{-1}\left]\right\}\times \\ & & \exp [-\frac{1}{2}\text{t}r\left\{\left(\right(B_{\xi }-{\tilde{B}}_{\xi }{)}^{\top }{K}_{\xi }(B_{\xi }-{\tilde{B}}_{\xi })){\Sigma }^{-1}\right\}]d{B}_{\xi }d\Sigma \pi \left(\xi \right)\\ & =& \left(2\pi {)}^{-\frac{|\xi |q}{2}}{|H_{\xi }|}^{-\frac{q}{2}}\right(2\pi {)}^{\frac{|\xi |q}{2}}|K_{\xi }{|}^{-\frac{q}{2}}\times \\ & & \int |\Sigma {|}^{\frac{|\xi |}{2}}|\Sigma {|}^{-\frac{1}{2}(n+|\xi |+\delta +q+1)}\exp \{-\frac{1}{2}\text{t}r[(Y^{\top }Y-{\tilde{B}}_{\xi }^{\top }{K}_{\xi }{\tilde{B}}_{\xi }+Q){\Sigma }^{-1}\left]\right\}d\Sigma \pi \left(\xi \right)\\ & =& {|H_{\xi }|}^{-\frac{q}{2}}|K_{\xi }{|}^{-\frac{q}{2}}\int |\Sigma {|}^{-\frac{1}{2}(n+\delta +q+1)}\times \\ & & \exp \{-\frac{1}{2}\text{t}r[(Y^{\top }Y-Y^{\top }{X}_{\xi }{K}_{\xi }^{-1}{X}_{\xi }^{\top }Y+Q){\Sigma }^{-1}\left]\right\}d\Sigma \pi \left(\xi \right)\\ & =& {|H_{\xi }{K}_{\xi }|}^{-\frac{q}{2}}\int |\Sigma {|}^{-\frac{1}{2}(n+\delta +q+1)}\times \\ & & \exp \{-\frac{1}{2}\text{t}r[\left(Y^{\top }\right(I_n-X_{\xi }{K}_{\xi }^{-1}{X}_{\xi }^{\top })Y+Q){\Sigma }^{-1}\left]\right\}d\Sigma \pi \left(\xi \right)\\ & \propto & {|H_{\xi }{K}_{\xi }|}^{-\frac{q}{2}}|Y^{\top }(I_n-X_{\xi }{K}_{\xi }^{-1}{X}_{\xi }^{\top })Y+Q{|}^{-\frac{n+\delta }{2}}\pi \left(\xi \right)\\ & \equiv & g\left(\xi |Y\right),& \text{(1)}\end{array}$$ where $K_{\xi }=X_{\xi }^{\top }{X}_{\xi }+H_{\xi }^{-1}$ and $g\left(\xi |Y\right)$ is the proportional form of $\pi \left(\xi |Y\right)$. (Brown and Vannucci 1998) set $Q=k{I}_q$ and $H_{\xi }=c(X_{\xi }^{\top }{X}_{\xi }{)}^{-1}$ and ${\omega }_j=\omega$, then $\pi \left(\xi \right)={\prod }_{j=1}^p{\omega }_j^{{\xi }_j}(1-{\omega }_j{)}^{1-{\xi }_j}={\omega }^{|\xi |}(1-\omega {)}^{p-|\xi |}$. Apply $Q=k{I}_q$ and $H_{\xi }=c(X_{\xi }^{\top }{X}_{\xi }{)}^{-1}$ into (1), hence we have $$\begin{array}{ccc}\log g\left(\xi |Y\right)& =& -\frac{n+\delta }{2}\log |k{I}_q+Y^{\top }Y-\frac{c}{c+1}{Y}^{\top }{X}_{\xi }\left(X_{\xi }^{\top }{X}_{\xi }^{-1}\right)X_{\xi }^{\top }Y|-\\ & & \frac{q|\xi |}{2}\log (c+1)+|\xi |\log \omega +(p-|\xi |)\log (1-\omega ).\end{array}$$ Exhaustive computation of $g\left(\xi |Y\right)$ for $2^p$ values of $\xi$ becomes prohibitive even in a super computer when $p$ is larger than 40. In this circumstances, using MCMC sampling to explore marginal posterior is possible. (Brown and Vannucci 1998) use Gibbs sampler to generate each $\xi$ value component-wise from full conditional distributions $\pi ({\xi }_j|{\xi }_{-j},Y)$, where ${\xi }_{-j}=\xi \setminus {\xi }_j$ . It is easy to show that $$\begin{array}{c}\pi ({\xi }_j=1|{\xi }_{-j},Y)=\frac{{\theta }_j}{{\theta }_j+1},\end{array}$$ where $$\begin{array}{c}{\theta }_j=\frac{g({\xi }_j=1,{\xi }_{-j}|Y)}{g({\xi }_j=0,{\xi }_{-j}|Y)}.\end{array}$$ The complete algorithm of (Brown and Vannucci 1998) is described in Algorithm 3:

For the hyperparameter, we follow the same setting with (Brown and Vannucci 1998), which are $\delta =2+q,\omega =20/p,k=2$ and $c=10$. Note that since our definition of inverse-wishart distribution is different from that of (Brown and Vannucci 1998), the values of $q$ and $k$ we set here are after adjustment to make them consistent with (Brown and Vannucci 1998).