Gibbs Linear Regression The variable selection problem has a natural Bayesian solution: Any collection of models having different sets of regressors can be computed via their Bayes factors.
A semiconjugate prior distribution
Let β ∼ M N ( β 0 , Σ 0 ) \beta \sim MN(\boldsymbol\beta_0, \Sigma_0) β ∼ MN ( β 0 , Σ 0 )
β ∣ y , X , σ 2 ∼ N ( ( Σ 0 − 1 + X T X / σ 2 ) − 1 ( Σ 0 − 1 β 0 + X T y / σ 2 ) , ( Σ 0 − 1 + X T X / σ 2 ) − 1 ) \beta\mid \mathbf y,\mathbf X, \sigma^2 \sim N((\Sigma_0^{-1}+\mathbf X^T\mathbf X/\sigma^2)^{-1}(\Sigma_0^{-1}\boldsymbol\beta_0+\mathbf X^T\mathbf y/\sigma^2), (\Sigma_0^{-1}+\mathbf X^T\mathbf X/\sigma^2)^{-1}) β ∣ y , X , σ 2 ∼ N (( Σ 0 − 1 + X T X / σ 2 ) − 1 ( Σ 0 − 1 β 0 + X T y / σ 2 ) , ( Σ 0 − 1 + X T X / σ 2 ) − 1 ) and let 1 / σ 2 ∼ G a ( ν 0 / 2 , ν 0 σ 0 2 / 2 ) 1/\sigma^2\sim Ga(\nu_0/2,\nu_0\sigma_0^2/2) 1/ σ 2 ∼ G a ( ν 0 /2 , ν 0 σ 0 2 /2 )
1 / σ 2 ∣ y , X , β ∼ G a ( [ ν 0 + n ] / 2 , [ ν 0 σ 0 2 + S S R ( β ) ] / 2 ) 1/\sigma^2\mid \mathbf y,\mathbf X,\boldsymbol\beta \sim Ga([\nu_0+n]/2, [\nu_0\sigma_0^2+SSR(\boldsymbol \beta)]/2) 1/ σ 2 ∣ y , X , β ∼ G a ([ ν 0 + n ] /2 , [ ν 0 σ 0 2 + SSR ( β )] /2 ) Then we can construct the following Gibbs sampler:
the parameter estimation should be invariant to changes in the scale of the regressors.
For the second case, we can derive a Monte Carlo approximation:
since
p ( σ 2 , β ∣ y , X ) ∝ p ( β ∣ y , X , σ 2 ) p ( σ 2 ∣ y , X ) p(\sigma^2,\boldsymbol\beta\mid \mathbf y, \mathbf X)\propto p(\boldsymbol\beta\mid \mathbf y,\mathbf X,\sigma^2)p(\sigma^2\mid \mathbf y,\mathbf X) p ( σ 2 , β ∣ y , X ) ∝ p ( β ∣ y , X , σ 2 ) p ( σ 2 ∣ y , X ) 
