From ad46fcb8cb6a36861682b674c59475f62b726476 Mon Sep 17 00:00:00 2001 From: Morten Hjorth-Jensen Date: Tue, 7 Nov 2017 10:52:06 +0100 Subject: [PATCH] chang to bayesian file --- doc/src/Bayesian/Bayesian.do.txt | 185 +++++++++++++++++++++++++++++++ 1 file changed, 185 insertions(+) diff --git a/doc/src/Bayesian/Bayesian.do.txt b/doc/src/Bayesian/Bayesian.do.txt index 06df75d36..be3a6fb3b 100644 --- a/doc/src/Bayesian/Bayesian.do.txt +++ b/doc/src/Bayesian/Bayesian.do.txt @@ -9,3 +9,188 @@ DATE: today !eblock +mathbf {y} -\mathbf {X} {\boldsymbol {\beta }})\right).} \rho(\mathbf{y}|\mathbf{X},\boldsymbol\beta,\sigma^{2}) \propto (\sigma^{2})^{-n/2} \exp\left(-\frac{1}{2{\sigma}^{2}}(\mathbf{y}- \mathbf{X} \boldsymbol\beta)^{\rm T}(\mathbf{y}- \mathbf{X} \boldsymbol\beta)\right).~ +The ordinary least squares solution is to estimate the coefficient vector using the Moore-Penrose pseudoinverse: + +{\displaystyle {\hat {\boldsymbol {\beta }}}=(\mathbf {X} ^{\rm {T}}\mathbf {X} )^{-1}\mathbf {X} ^{\rm {T}}\mathbf {y} } \hat{\boldsymbol\beta} = (\mathbf{X}^{\rm T}\mathbf{X})^{-1}\mathbf{X}^{\rm T}\mathbf{y} +where {\displaystyle \mathbf {X} } \mathbf {X} is the {\displaystyle n\times k} n \times k design matrix, each row of which is a predictor vector {\displaystyle \mathbf {x} _{i}^{\rm {T}}} \mathbf{x}_{i}^{\rm T}; and {\displaystyle \mathbf {y} } \mathbf {y} is the column {\displaystyle n} n-vector {\displaystyle [y_{1}\;\cdots \;y_{n}]^{\rm {T}}} [y_1 \; \cdots \; y_n]^{\rm T}. + +This is a frequentist approach, and it assumes that there are enough measurements to say something meaningful about {\displaystyle {\boldsymbol {\beta }}} {\boldsymbol {\beta }}. In the Bayesian approach, the data are supplemented with additional information in the form of a prior probability distribution. The prior belief about the parameters is combined with the data's likelihood function according to Bayes theorem to yield the posterior belief about the parameters {\displaystyle {\boldsymbol {\beta }}} {\boldsymbol {\beta }} and {\displaystyle \sigma } \sigma . The prior can take different functional forms depending on the domain and the information that is available in a priori. + +With conjugate priors[edit] +Conjugate prior distribution[edit] +For an arbitrary prior distribution, there may be no analytical solution for the posterior distribution. In this section, we will consider a so-called conjugate prior for which the posterior distribution can be derived analytically. + +A prior {\displaystyle \rho ({\boldsymbol {\beta }},\sigma ^{2})} \rho(\boldsymbol\beta,\sigma^{2}) is conjugate to this likelihood function if it has the same functional form with respect to {\displaystyle {\boldsymbol {\beta }}} {\boldsymbol {\beta }} and {\displaystyle \sigma } \sigma . Since the log-likelihood is quadratic in {\displaystyle {\boldsymbol {\beta }}} {\boldsymbol {\beta }}, the log-likelihood is re-written such that the likelihood becomes normal in {\displaystyle ({\boldsymbol {\beta }}-{\hat {\boldsymbol {\beta }}})} (\boldsymbol\beta-\hat{\boldsymbol\beta}). Write + +{\displaystyle (\mathbf {y} -\mathbf {X} {\boldsymbol {\beta }})^{\rm {T}}(\mathbf {y} -\mathbf {X} {\boldsymbol {\beta }})=(\mathbf {y} -\mathbf {X} {\hat {\boldsymbol {\beta }}})^{\rm {T}}(\mathbf {y} -\mathbf {X} {\hat {\boldsymbol {\beta }}})+({\boldsymbol {\beta }}-{\hat {\boldsymbol {\beta }}})^{\rm {T}}(\mathbf {X} ^{\rm {T}}\mathbf {X} )({\boldsymbol {\beta }}-{\hat {\boldsymbol {\beta }}}).} {\displaystyle (\mathbf {y} -\mathbf {X} {\boldsymbol {\beta }})^{\rm {T}}(\mathbf {y} -\mathbf {X} {\boldsymbol {\beta }})=(\mathbf {y} -\mathbf {X} {\hat {\boldsymbol {\beta }}})^{\rm {T}}(\mathbf {y} -\mathbf {X} {\hat {\boldsymbol {\beta }}})+({\boldsymbol {\beta }}-{\hat {\boldsymbol {\beta }}})^{\rm {T}}(\mathbf {X} ^{\rm {T}}\mathbf {X} )({\boldsymbol {\beta }}-{\hat {\boldsymbol {\beta }}}).} +The likelihood is now re-written as + +{\displaystyle \rho (\mathbf {y} |\mathbf {X} ,{\boldsymbol {\beta }},\sigma ^{2})\propto (\sigma ^{2})^{-{\frac {v}{2}}}\exp \left(-{\frac {vs^{2}}{2{\sigma }^{2}}}\right)(\sigma ^{2})^{-{\frac {n-v}{2}}}\exp \left(-{\frac {1}{2{\sigma }^{2}}}({\boldsymbol {\beta }}-{\hat {\boldsymbol {\beta }}})^{\rm {T}}(\mathbf {X} ^{\rm {T}}\mathbf {X} )({\boldsymbol {\beta }}-{\hat {\boldsymbol {\beta }}})\right),} {\displaystyle \rho (\mathbf {y} |\mathbf {X} ,{\boldsymbol {\beta }},\sigma ^{2})\propto (\sigma ^{2})^{-{\frac {v}{2}}}\exp \left(-{\frac {vs^{2}}{2{\sigma }^{2}}}\right)(\sigma ^{2})^{-{\frac {n-v}{2}}}\exp \left(-{\frac {1}{2{\sigma }^{2}}}({\boldsymbol {\beta }}-{\hat {\boldsymbol {\beta }}})^{\rm {T}}(\mathbf {X} ^{\rm {T}}\mathbf {X} )({\boldsymbol {\beta }}-{\hat {\boldsymbol {\beta }}})\right),} +where + +{\displaystyle vs^{2}=(\mathbf {y} -\mathbf {X} {\hat {\boldsymbol {\beta }}})^{\rm {T}}(\mathbf {y} -\mathbf {X} {\hat {\boldsymbol {\beta }}})\quad {\text{ and }}\quad v=n-k,} {\displaystyle vs^{2}=(\mathbf {y} -\mathbf {X} {\hat {\boldsymbol {\beta }}})^{\rm {T}}(\mathbf {y} -\mathbf {X} {\hat {\boldsymbol {\beta }}})\quad {\text{ and }}\quad v=n-k,} +where {\displaystyle k} k is the number of regression coefficients. + +This suggests a form for the prior: + +{\displaystyle \rho ({\boldsymbol {\beta }},\sigma ^{2})=\rho (\sigma ^{2})\rho ({\boldsymbol {\beta }}|\sigma ^{2}),} \rho(\boldsymbol\beta,\sigma^{2}) = \rho(\sigma^{2})\rho(\boldsymbol\beta|\sigma^{2}), +where {\displaystyle \rho (\sigma ^{2})} \rho(\sigma^{2}) is an inverse-gamma distribution + +{\displaystyle \rho (\sigma ^{2})\propto (\sigma ^{2})^{-{\frac {v_{0}}{2}}-1}\exp \left(-{\frac {v_{0}s_{0}^{2}}{2{\sigma }^{2}}}\right).} {\displaystyle \rho (\sigma ^{2})\propto (\sigma ^{2})^{-{\frac {v_{0}}{2}}-1}\exp \left(-{\frac {v_{0}s_{0}^{2}}{2{\sigma }^{2}}}\right).} +In the notation introduced in the inverse-gamma distribution article, this is the density of an {\displaystyle {\text{Inv-Gamma}}(a_{0},b_{0})} \text{Inv-Gamma}( a_0,b_0) distribution with {\displaystyle a_{0}={\tfrac {v_{0}}{2}}} {\displaystyle a_{0}={\tfrac {v_{0}}{2}}} and {\displaystyle b_{0}={\tfrac {1}{2}}v_{0}s_{0}^{2}} {\displaystyle b_{0}={\tfrac {1}{2}}v_{0}s_{0}^{2}} with {\displaystyle v_{0}} v_{0} and {\displaystyle s_{0}^{2}} s_{0}^{2} as the prior values of {\displaystyle v} v and {\displaystyle s^{2}} s^{2}, respectively. Equivalently, it can also be described as a scaled inverse chi-squared distribution, {\displaystyle {\mbox{Scale-inv-}}\chi ^{2}(v_{0},s_{0}^{2}).} \mbox{Scale-inv-}\chi^2(v_0, s_0^2). + +Further the conditional prior density {\displaystyle \rho ({\boldsymbol {\beta }}|\sigma ^{2})} \rho(\boldsymbol\beta|\sigma^{2}) is a normal distribution, + +{\displaystyle \rho ({\boldsymbol {\beta }}|\sigma ^{2})\propto (\sigma ^{2})^{-{\frac {k}{2}}}\exp \left(-{\frac {1}{2{\sigma }^{2}}}({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{0})^{\rm {T}}\mathbf {\Lambda } _{0}({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{0})\right).} {\displaystyle \rho ({\boldsymbol {\beta }}|\sigma ^{2})\propto (\sigma ^{2})^{-{\frac {k}{2}}}\exp \left(-{\frac {1}{2{\sigma }^{2}}}({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{0})^{\rm {T}}\mathbf {\Lambda } _{0}({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{0})\right).} +In the notation of the normal distribution, the conditional prior distribution is {\displaystyle {\mathcal {N}}\left({\boldsymbol {\mu }}_{0},\sigma ^{2}\mathbf {\Lambda } _{0}^{-1}\right).} \mathcal{N}\left(\boldsymbol\mu_0, \sigma^2\mathbf{\Lambda}_0^{-1}\right). + +Posterior distribution[edit] +With the prior now specified, the posterior distribution can be expressed as + +{\displaystyle {\begin{aligned}\rho ({\boldsymbol {\beta }},\sigma ^{2}|\mathbf {y} ,\mathbf {X} )&\propto \rho (\mathbf {y} |\mathbf {X} ,{\boldsymbol {\beta }},\sigma ^{2})\rho ({\boldsymbol {\beta }}|\sigma ^{2})\rho (\sigma ^{2})\\&\propto (\sigma ^{2})^{-{\frac {n}{2}}}\exp \left(-{\frac {1}{2{\sigma }^{2}}}(\mathbf {y} -\mathbf {X} {\boldsymbol {\beta }})^{\rm {T}}(\mathbf {y} -\mathbf {X} {\boldsymbol {\beta }})\right)(\sigma ^{2})^{-{\frac {k}{2}}}\exp \left(-{\frac {1}{2{\sigma }^{2}}}({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{0})^{\rm {T}}{\boldsymbol {\Lambda }}_{0}({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{0})\right)(\sigma ^{2})^{-(a_{0}+1)}\exp \left(-{\frac {b_{0}}{{\sigma }^{2}}}\right)\end{aligned}}} {\displaystyle {\begin{aligned}\rho ({\boldsymbol {\beta }},\sigma ^{2}|\mathbf {y} ,\mathbf {X} )&\propto \rho (\mathbf {y} |\mathbf {X} ,{\boldsymbol {\beta }},\sigma ^{2})\rho ({\boldsymbol {\beta }}|\sigma ^{2})\rho (\sigma ^{2})\\&\propto (\sigma ^{2})^{-{\frac {n}{2}}}\exp \left(-{\frac {1}{2{\sigma }^{2}}}(\mathbf {y} -\mathbf {X} {\boldsymbol {\beta }})^{\rm {T}}(\mathbf {y} -\mathbf {X} {\boldsymbol {\beta }})\right)(\sigma ^{2})^{-{\frac {k}{2}}}\exp \left(-{\frac {1}{2{\sigma }^{2}}}({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{0})^{\rm {T}}{\boldsymbol {\Lambda }}_{0}({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{0})\right)(\sigma ^{2})^{-(a_{0}+1)}\exp \left(-{\frac {b_{0}}{{\sigma }^{2}}}\right)\end{aligned}}} +With some re-arrangement,[1] the posterior can be re-written so that the posterior mean {\displaystyle {\boldsymbol {\mu }}_{n}} \boldsymbol\mu_n of the parameter vector {\displaystyle {\boldsymbol {\beta }}} {\boldsymbol {\beta }} can be expressed in terms of the least squares estimator {\displaystyle {\hat {\boldsymbol {\beta }}}} \hat{\boldsymbol\beta} and the prior mean {\displaystyle {\boldsymbol {\mu }}_{0}} {\boldsymbol {\mu }}_{0}, with the strength of the prior indicated by the prior precision matrix {\displaystyle {\boldsymbol {\Lambda }}_{0}} \boldsymbol\Lambda_0 + +{\displaystyle {\boldsymbol {\mu }}_{n}=(\mathbf {X} ^{\rm {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0})^{-1}(\mathbf {X} ^{\rm {T}}\mathbf {X} {\hat {\boldsymbol {\beta }}}+{\boldsymbol {\Lambda }}_{0}{\boldsymbol {\mu }}_{0}).} \boldsymbol\mu_n = (\mathbf{X}^{\rm T}\mathbf{X}+\boldsymbol\Lambda_0)^{-1}(\mathbf{X}^{\rm T}\mathbf{X}\hat{\boldsymbol\beta}+\boldsymbol\Lambda_0\boldsymbol\mu_0) . +To justify that {\displaystyle {\boldsymbol {\mu }}_{n}} \boldsymbol\mu_n is indeed the posterior mean, the quadratic terms in the exponential can be re-arranged as a quadratic form in {\displaystyle {\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{n}} \boldsymbol\beta-\boldsymbol\mu_n.[2] + +{\displaystyle (\mathbf {y} -\mathbf {X} {\boldsymbol {\beta }})^{\rm {T}}(\mathbf {y} -\mathbf {X} {\boldsymbol {\beta }})+({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{0})^{\rm {T}}{\boldsymbol {\Lambda }}_{0}({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{0})=({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{n})^{\rm {T}}(\mathbf {X} ^{\rm {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0})({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{n})+\mathbf {y} ^{\rm {T}}\mathbf {y} -{\boldsymbol {\mu }}_{n}^{\rm {T}}(\mathbf {X} ^{\rm {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0}){\boldsymbol {\mu }}_{n}+{\boldsymbol {\mu }}_{0}^{\rm {T}}{\boldsymbol {\Lambda }}_{0}{\boldsymbol {\mu }}_{0}.} {\displaystyle (\mathbf {y} -\mathbf {X} {\boldsymbol {\beta }})^{\rm {T}}(\mathbf {y} -\mathbf {X} {\boldsymbol {\beta }})+({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{0})^{\rm {T}}{\boldsymbol {\Lambda }}_{0}({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{0})=({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{n})^{\rm {T}}(\mathbf {X} ^{\rm {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0})({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{n})+\mathbf {y} ^{\rm {T}}\mathbf {y} -{\boldsymbol {\mu }}_{n}^{\rm {T}}(\mathbf {X} ^{\rm {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0}){\boldsymbol {\mu }}_{n}+{\boldsymbol {\mu }}_{0}^{\rm {T}}{\boldsymbol {\Lambda }}_{0}{\boldsymbol {\mu }}_{0}.} +Now the posterior can be expressed as a normal distribution times an inverse-gamma distribution: + +{\displaystyle \rho ({\boldsymbol {\beta }},\sigma ^{2}|\mathbf {y} ,\mathbf {X} )\propto (\sigma ^{2})^{-{\frac {k}{2}}}\exp \left(-{\frac {1}{2{\sigma }^{2}}}({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{n})^{\rm {T}}(\mathbf {X} ^{\rm {T}}\mathbf {X} +\mathbf {\Lambda } _{0})({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{n})\right)(\sigma ^{2})^{-{\frac {n+2a_{0}}{2}}-1}\exp \left(-{\frac {2b_{0}+\mathbf {y} ^{\rm {T}}\mathbf {y} -{\boldsymbol {\mu }}_{n}^{\rm {T}}(\mathbf {X} ^{\rm {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0}){\boldsymbol {\mu }}_{n}+{\boldsymbol {\mu }}_{0}^{\rm {T}}{\boldsymbol {\Lambda }}_{0}{\boldsymbol {\mu }}_{0}}{2{\sigma }^{2}}}\right).} {\displaystyle \rho ({\boldsymbol {\beta }},\sigma ^{2}|\mathbf {y} ,\mathbf {X} )\propto (\sigma ^{2})^{-{\frac {k}{2}}}\exp \left(-{\frac {1}{2{\sigma }^{2}}}({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{n})^{\rm {T}}(\mathbf {X} ^{\rm {T}}\mathbf {X} +\mathbf {\Lambda } _{0})({\boldsymbol {\beta }}-{\boldsymbol {\mu }}_{n})\right)(\sigma ^{2})^{-{\frac {n+2a_{0}}{2}}-1}\exp \left(-{\frac {2b_{0}+\mathbf {y} ^{\rm {T}}\mathbf {y} -{\boldsymbol {\mu }}_{n}^{\rm {T}}(\mathbf {X} ^{\rm {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0}){\boldsymbol {\mu }}_{n}+{\boldsymbol {\mu }}_{0}^{\rm {T}}{\boldsymbol {\Lambda }}_{0}{\boldsymbol {\mu }}_{0}}{2{\sigma }^{2}}}\right).} +Therefore, the posterior distribution can be parametrized as follows. + +{\displaystyle \rho ({\boldsymbol {\beta }},\sigma ^{2}|\mathbf {y} ,\mathbf {X} )\propto \rho ({\boldsymbol {\beta }}|\sigma ^{2},\mathbf {y} ,\mathbf {X} )\rho (\sigma ^{2}|\mathbf {y} ,\mathbf {X} ),} \rho(\boldsymbol\beta,\sigma^{2}|\mathbf{y},\mathbf{X}) \propto \rho(\boldsymbol\beta |\sigma^{2},\mathbf{y},\mathbf{X}) \rho(\sigma^{2}|\mathbf{y},\mathbf{X}), +where the two factors correspond to the densities of {\displaystyle {\mathcal {N}}\left({\boldsymbol {\mu }}_{n},\sigma ^{2}{\boldsymbol {\Lambda }}_{n}^{-1}\right)\,} \mathcal{N}\left( \boldsymbol\mu_n, \sigma^2\boldsymbol\Lambda_n^{-1} \right)\, and {\displaystyle {\text{Inv-Gamma}}\left(a_{n},b_{n}\right)} \text{Inv-Gamma}\left(a_n,b_n \right) distributions, with the parameters of these given by + +{\displaystyle {\boldsymbol {\Lambda }}_{n}=(\mathbf {X} ^{\rm {T}}\mathbf {X} +\mathbf {\Lambda } _{0}),\quad {\boldsymbol {\mu }}_{n}=({\boldsymbol {\Lambda }}_{n})^{-1}(\mathbf {X} ^{\rm {T}}\mathbf {X} {\hat {\boldsymbol {\beta }}}+{\boldsymbol {\Lambda }}_{0}{\boldsymbol {\mu }}_{0}),} \boldsymbol\Lambda_n=(\mathbf{X}^{\rm T}\mathbf{X}+\mathbf{\Lambda}_0), \quad \boldsymbol\mu_n = (\boldsymbol\Lambda_n)^{-1}(\mathbf{X}^{\rm T}\mathbf{X}\hat{\boldsymbol\beta}+\boldsymbol\Lambda_0\boldsymbol\mu_0) , +{\displaystyle a_{n}=a_{0}+{\frac {n}{2}},\qquad b_{n}=b_{0}+{\frac {1}{2}}(\mathbf {y} ^{\rm {T}}\mathbf {y} +{\boldsymbol {\mu }}_{0}^{\rm {T}}{\boldsymbol {\Lambda }}_{0}{\boldsymbol {\mu }}_{0}-{\boldsymbol {\mu }}_{n}^{\rm {T}}{\boldsymbol {\Lambda }}_{n}{\boldsymbol {\mu }}_{n}).} a_{n}=a_{0}+{\frac {n}{2}},\qquad b_{n}=b_{0}+{\frac {1}{2}}({\mathbf {y}}^{{{\rm {T}}}}{\mathbf {y}}+{\boldsymbol \mu }_{0}^{{{\rm {T}}}}{\boldsymbol \Lambda }_{0}{\boldsymbol \mu }_{0}-{\boldsymbol \mu }_{n}^{{{\rm {T}}}}{\boldsymbol \Lambda }_{n}{\boldsymbol \mu }_{n}). +This can be interpreted as Bayesian learning where the parameters are updated according to the following equations. + +{\displaystyle {\boldsymbol {\mu }}_{n}=(\mathbf {X} ^{\rm {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0})^{-1}({\boldsymbol {\Lambda }}_{0}{\boldsymbol {\mu }}_{0}+\mathbf {X} ^{\rm {T}}\mathbf {X} {\hat {\boldsymbol {\beta }}})=(\mathbf {X} ^{\rm {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0})^{-1}({\boldsymbol {\Lambda }}_{0}{\boldsymbol {\mu }}_{0}+\mathbf {X} ^{\rm {T}}\mathbf {y} ),} \boldsymbol\mu_n=(\mathbf{X}^{\rm T}\mathbf{X}+\boldsymbol\Lambda_0)^{-1} (\boldsymbol\Lambda_0\boldsymbol\mu_0+\mathbf{X}^{\rm T}\mathbf{X}\hat{\boldsymbol\beta})=(\mathbf{X}^{\rm T}\mathbf{X}+\boldsymbol\Lambda_0)^{-1} (\boldsymbol\Lambda_0\boldsymbol\mu_0+\mathbf{X}^{\rm T}\mathbf{y}) , +{\displaystyle {\boldsymbol {\Lambda }}_{n}=(\mathbf {X} ^{\rm {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0}),} \boldsymbol\Lambda_n=(\mathbf{X}^{\rm T}\mathbf{X}+\boldsymbol\Lambda_0) , +{\displaystyle a_{n}=a_{0}+{\frac {n}{2}},} a_{n}=a_{0}+{\frac {n}{2}}, +{\displaystyle b_{n}=b_{0}+{\frac {1}{2}}(\mathbf {y} ^{\rm {T}}\mathbf {y} +{\boldsymbol {\mu }}_{0}^{\rm {T}}{\boldsymbol {\Lambda }}_{0}{\boldsymbol {\mu }}_{0}-{\boldsymbol {\mu }}_{n}^{\rm {T}}{\boldsymbol {\Lambda }}_{n}{\boldsymbol {\mu }}_{n}).} b_{n}=b_{0}+{\frac {1}{2}}({\mathbf {y}}^{{{\rm {T}}}}{\mathbf {y}}+{\boldsymbol \mu }_{0}^{{{\rm {T}}}}{\boldsymbol \Lambda }_{0}{\boldsymbol \mu }_{0}-{\boldsymbol \mu }_{n}^{{{\rm {T}}}}{\boldsymbol \Lambda }_{n}{\boldsymbol \mu }_{n}). +Model evidence[edit] +The model evidence {\displaystyle p(\mathbf {y} |m)} p(\mathbf{y}|m) is the probability of the data given the model {\displaystyle m} m. It is also known as the marginal likelihood, and as the prior predictive density. Here, the model is defined by the likelihood function {\displaystyle p(\mathbf {y} |\mathbf {X} ,{\boldsymbol {\beta }},\sigma )} p(\mathbf{y}|\mathbf{X},\boldsymbol\beta,\sigma) and the prior distribution on the parameters, i.e. {\displaystyle p({\boldsymbol {\beta }},\sigma )} p(\boldsymbol\beta,\sigma). The model evidence captures in a single number how well such a model explains the observations. The model evidence of the Bayesian linear regression model presented in this section can be used to compare competing linear models by Bayesian model comparison. These models may differ in the number and values of the predictor variables as well as in their priors on the model parameters. Model complexity is already taken into account by the model evidence, because it marginalizes out the parameters by integrating {\displaystyle p(\mathbf {y} ,{\boldsymbol {\beta }},\sigma |\mathbf {X} )} p(\mathbf{y},\boldsymbol\beta,\sigma|\mathbf{X}) over all possible values of {\displaystyle {\boldsymbol {\beta }}} {\boldsymbol {\beta }} and {\displaystyle \sigma } \sigma . + +{\displaystyle p(\mathbf {y} |m)=\int p(\mathbf {y} |\mathbf {X} ,{\boldsymbol {\beta }},\sigma )\,p({\boldsymbol {\beta }},\sigma )\,d{\boldsymbol {\beta }}\,d\sigma } p(\mathbf{y}|m)=\int p(\mathbf{y}|\mathbf{X},\boldsymbol\beta,\sigma)\, p(\boldsymbol\beta,\sigma)\, d\boldsymbol\beta\, d\sigma +This integral can be computed analytically and the solution is given in the following equation.[3] + +{\displaystyle p(\mathbf {y} |m)={\frac {1}{(2\pi )^{n/2}}}{\sqrt {\frac {\det({\boldsymbol {\Lambda }}_{0})}{\det({\boldsymbol {\Lambda }}_{n})}}}\cdot {\frac {b_{0}^{a_{0}}}{b_{n}^{a_{n}}}}\cdot {\frac {\Gamma (a_{n})}{\Gamma (a_{0})}}} p(\mathbf{y}|m)=\frac{1}{(2\pi)^{n/2}}\sqrt{\frac{\det(\boldsymbol\Lambda_0)}{\det(\boldsymbol\Lambda_n)}} \cdot \frac{b_0^{a_0}}{b_n^{a_n}} \cdot \frac{\Gamma(a_n)}{\Gamma(a_0)} +Here {\displaystyle \Gamma } \Gamma denotes the gamma function. Because we have chosen a conjugate prior, the marginal likelihood can also be easily computed by evaluating the following equality for arbitrary values of {\displaystyle {\boldsymbol {\beta }}} {\boldsymbol {\beta }} and {\displaystyle \sigma } \sigma . + +{\displaystyle p(\mathbf {y} |m)={\frac {p({\boldsymbol {\beta }},\sigma |m)\,p(\mathbf {y} |\mathbf {X} ,{\boldsymbol {\beta }},\sigma ,m)}{p({\boldsymbol {\beta }},\sigma |\mathbf {y} ,\mathbf {X} ,m)}}} p(\mathbf{y}|m)=\frac{p(\boldsymbol\beta,\sigma|m)\, p(\mathbf{y}|\mathbf{X},\boldsymbol\beta,\sigma,m)}{p(\boldsymbol\beta,\sigma|\mathbf{y},\mathbf{X},m)} +Note that this equation is nothing but a re-arrangement of Bayes theorem. Inserting the formulas for the prior, the likelihood, and the posterior and simplifying the resulting expression leads to the analytic expression given above. + +Other cases[edit] +In general, it may be impossible or impractical to derive the posterior distribution analytically. However, it is possible to approximate the posterior by an approximate Bayesian inference method such as Monte Carlo sampling[4] or variational Bayes. + +The special case {\displaystyle {\boldsymbol {\mu }}_{0}=0,\mathbf {\Lambda } _{0}=c\mathbf {I} } \boldsymbol\mu_0=0, \mathbf{\Lambda}_0 = c\mathbf{I} is called ridge regression. + +A similar analysis can be performed for the general case of the multivariate regression and part of this provides for Bayesian estimation of covariance matrices: see Bayesian multivariate linear regression.~ + + + + + + + + +ts ,\epsilon _{i,m}\}} \{\epsilon _{{i,1}},\ldots ,\epsilon _{{i,m}}\} are all correlated. Equivalently, it can be viewed as a single regression problem where the outcome is a row vector {\displaystyle \mathbf {y} _{i}^{\rm {T}}} {\mathbf {y}}_{i}^{{{\rm {T}}}} and the regression coefficient vectors are stacked next to each other, as follows: + +{\displaystyle \mathbf {y} _{i}^{\rm {T}}=\mathbf {x} _{i}^{\rm {T}}\mathbf {B} +{\boldsymbol {\epsilon }}_{i}^{\rm {T}}.} {\mathbf {y}}_{i}^{{{\rm {T}}}}={\mathbf {x}}_{i}^{{{\rm {T}}}}{\mathbf {B}}+{\boldsymbol \epsilon }_{{i}}^{{{\rm {T}}}}. +The coefficient matrix B is a {\displaystyle k\times m} k\times m matrix where the coefficient vectors {\displaystyle {\boldsymbol {\beta }}_{1},\ldots ,{\boldsymbol {\beta }}_{m}} {\boldsymbol \beta }_{1},\ldots ,{\boldsymbol \beta }_{m} for each regression problem are stacked horizontally: + +{\displaystyle \mathbf {B} ={\begin{bmatrix}{\begin{pmatrix}\\{\boldsymbol {\beta }}_{1}\\\\\end{pmatrix}}\cdots {\begin{pmatrix}\\{\boldsymbol {\beta }}_{m}\\\\\end{pmatrix}}\end{bmatrix}}={\begin{bmatrix}{\begin{pmatrix}\beta _{1,1}\\\vdots \\\beta _{k,1}\\\end{pmatrix}}\cdots {\begin{pmatrix}\beta _{1,m}\\\vdots \\\beta _{k,m}\\\end{pmatrix}}\end{bmatrix}}.} {\mathbf {B}}={\begin{bmatrix}{\begin{pmatrix}\\{\boldsymbol \beta }_{1}\\\\\end{pmatrix}}\cdots {\begin{pmatrix}\\{\boldsymbol \beta }_{m}\\\\\end{pmatrix}}\end{bmatrix}}={\begin{bmatrix}{\begin{pmatrix}\beta _{{1,1}}\\\vdots \\\beta _{{k,1}}\\\end{pmatrix}}\cdots {\begin{pmatrix}\beta _{{1,m}}\\\vdots \\\beta _{{k,m}}\\\end{pmatrix}}\end{bmatrix}}. +The noise vector {\displaystyle {\boldsymbol {\epsilon }}_{i}} {\boldsymbol \epsilon }_{{i}} for each observation i is jointly normal, so that the outcomes for a given observation are correlated: + +{\displaystyle {\boldsymbol {\epsilon }}_{i}\sim N(0,{\boldsymbol {\Sigma }}_{\epsilon }).} {\displaystyle {\boldsymbol {\epsilon }}_{i}\sim N(0,{\boldsymbol {\Sigma }}_{\epsilon }).} +We can write the entire regression problem in matrix form as: + +{\displaystyle \mathbf {Y} =\mathbf {X} \mathbf {B} +\mathbf {E} ,} {\mathbf {Y}}={\mathbf {X}}{\mathbf {B}}+{\mathbf {E}}, +where Y and E are {\displaystyle n\times m} n\times m matrices. The design matrix X is an {\displaystyle n\times k} n \times k matrix with the observations stacked vertically, as in the standard linear regression setup: + +{\displaystyle \mathbf {X} ={\begin{bmatrix}\mathbf {x} _{1}^{\rm {T}}\\\mathbf {x} _{2}^{\rm {T}}\\\vdots \\\mathbf {x} _{n}^{\rm {T}}\end{bmatrix}}={\begin{bmatrix}x_{1,1}&\cdots &x_{1,k}\\x_{2,1}&\cdots &x_{2,k}\\\vdots &\ddots &\vdots \\x_{n,1}&\cdots &x_{n,k}\end{bmatrix}}.} {\mathbf {X}}={\begin{bmatrix}{\mathbf {x}}_{1}^{{{\rm {T}}}}\\{\mathbf {x}}_{2}^{{{\rm {T}}}}\\\vdots \\{\mathbf {x}}_{n}^{{{\rm {T}}}}\end{bmatrix}}={\begin{bmatrix}x_{{1,1}}&\cdots &x_{{1,k}}\\x_{{2,1}}&\cdots &x_{{2,k}}\\\vdots &\ddots &\vdots \\x_{{n,1}}&\cdots &x_{{n,k}}\end{bmatrix}}. +The classical, frequentists linear least squares solution is to simply estimate the matrix of regression coefficients {\displaystyle {\hat {\mathbf {B} }}} {\hat {{\mathbf {B}}}} using the Moore-Penrose pseudoinverse: + +{\displaystyle {\hat {\mathbf {B} }}=(\mathbf {X} ^{\rm {T}}\mathbf {X} )^{-1}\mathbf {X} ^{\rm {T}}\mathbf {Y} } {\hat {{\mathbf {B}}}}=({\mathbf {X}}^{{{\rm {T}}}}{\mathbf {X}})^{{-1}}{\mathbf {X}}^{{{\rm {T}}}}{\mathbf {Y}}. +To obtain the Bayesian solution, we need to specify the conditional likelihood and then find the appropriate conjugate prior. As with the univariate case of linear Bayesian regression, we will find that we can specify a natural conditional conjugate prior (which is scale dependent). + +Let us write our conditional likelihood as[1] + +{\displaystyle \rho (\mathbf {E} |{\boldsymbol {\Sigma }}_{\epsilon })\propto |{\boldsymbol {\Sigma }}_{\epsilon }|^{-n/2}\exp(-{\frac {1}{2}}{\rm {tr}}(\mathbf {E} ^{\rm {T}}\mathbf {E} {\boldsymbol {\Sigma }}_{\epsilon }^{-1})),} {\displaystyle \rho (\mathbf {E} |{\boldsymbol {\Sigma }}_{\epsilon })\propto |{\boldsymbol {\Sigma }}_{\epsilon }|^{-n/2}\exp(-{\frac {1}{2}}{\rm {tr}}(\mathbf {E} ^{\rm {T}}\mathbf {E} {\boldsymbol {\Sigma }}_{\epsilon }^{-1})),} +writing the error {\displaystyle \mathbf {E} } \mathbf {E} in terms of {\displaystyle \mathbf {Y} ,\mathbf {X} ,} {\mathbf {Y}},{\mathbf {X}}, and {\displaystyle \mathbf {B} } \mathbf {B} yields + +{\displaystyle \rho (\mathbf {Y} |\mathbf {X} ,\mathbf {B} ,{\boldsymbol {\Sigma }}_{\epsilon })\propto |{\boldsymbol {\Sigma }}_{\epsilon }|^{-n/2}\exp(-{\frac {1}{2}}{\rm {tr}}((\mathbf {Y} -\mathbf {X} \mathbf {\mathbf {B} } )^{\rm {T}}(\mathbf {Y} -\mathbf {X} \mathbf {\mathbf {B} } ){\boldsymbol {\Sigma }}_{\epsilon }^{-1})),} {\displaystyle \rho (\mathbf {Y} |\mathbf {X} ,\mathbf {B} ,{\boldsymbol {\Sigma }}_{\epsilon })\propto |{\boldsymbol {\Sigma }}_{\epsilon }|^{-n/2}\exp(-{\frac {1}{2}}{\rm {tr}}((\mathbf {Y} -\mathbf {X} \mathbf {\mathbf {B} } )^{\rm {T}}(\mathbf {Y} -\mathbf {X} \mathbf {\mathbf {B} } ){\boldsymbol {\Sigma }}_{\epsilon }^{-1})),} +We seek a natural conjugate prior—a joint density {\displaystyle \rho (\mathbf {B} ,\Sigma _{\epsilon })} \rho ({\mathbf {B}},\Sigma _{{\epsilon }}) which is of the same functional form as the likelihood. Since the likelihood is quadratic in {\displaystyle \mathbf {B} } \mathbf {B} , we re-write the likelihood so it is normal in {\displaystyle (\mathbf {B} -{\hat {\mathbf {B} }})} ({\mathbf {B}}-{\hat {{\mathbf {B}}}}) (the deviation from classical sample estimate). + +Using the same technique as with Bayesian linear regression, we decompose the exponential term using a matrix-form of the sum-of-squares technique. Here, however, we will also need to use the Matrix Differential Calculus (Kronecker product and vectorization transformations). + +First, let us apply sum-of-squares to obtain new expression for the likelihood: + +{\displaystyle \rho (\mathbf {Y} |\mathbf {X} ,\mathbf {B} ,{\boldsymbol {\Sigma }}_{\epsilon })\propto |{\boldsymbol {\Sigma }}_{\epsilon }|^{-(n-k)/2}\exp(-{\rm {tr}}({\frac {1}{2}}\mathbf {S} ^{\rm {T}}\mathbf {S} {\boldsymbol {\Sigma }}_{\epsilon }^{-1}))|{\boldsymbol {\Sigma }}_{\epsilon }|^{-k/2}\exp(-{\frac {1}{2}}{\rm {tr}}((\mathbf {B} -{\hat {\mathbf {B} }})^{\rm {T}}\mathbf {X} ^{\rm {T}}\mathbf {X} (\mathbf {B} -{\hat {\mathbf {B} }}){\boldsymbol {\Sigma }}_{\epsilon }^{-1})),} {\displaystyle \rho (\mathbf {Y} |\mathbf {X} ,\mathbf {B} ,{\boldsymbol {\Sigma }}_{\epsilon })\propto |{\boldsymbol {\Sigma }}_{\epsilon }|^{-(n-k)/2}\exp(-{\rm {tr}}({\frac {1}{2}}\mathbf {S} ^{\rm {T}}\mathbf {S} {\boldsymbol {\Sigma }}_{\epsilon }^{-1}))|{\boldsymbol {\Sigma }}_{\epsilon }|^{-k/2}\exp(-{\frac {1}{2}}{\rm {tr}}((\mathbf {B} -{\hat {\mathbf {B} }})^{\rm {T}}\mathbf {X} ^{\rm {T}}\mathbf {X} (\mathbf {B} -{\hat {\mathbf {B} }}){\boldsymbol {\Sigma }}_{\epsilon }^{-1})),} +{\displaystyle \mathbf {S} =\mathbf {Y} -\mathbf {X} {\hat {\mathbf {B} }}} {\displaystyle \mathbf {S} =\mathbf {Y} -\mathbf {X} {\hat {\mathbf {B} }}} +We would like to develop a conditional form for the priors: + +{\displaystyle \rho (\mathbf {B} ,{\boldsymbol {\Sigma }}_{\epsilon })=\rho ({\boldsymbol {\Sigma }}_{\epsilon })\rho (\mathbf {B} |{\boldsymbol {\Sigma }}_{\epsilon }),} \rho ({\mathbf {B}},{\boldsymbol \Sigma }_{{\epsilon }})=\rho ({\boldsymbol \Sigma }_{{\epsilon }})\rho ({\mathbf {B}}|{\boldsymbol \Sigma }_{{\epsilon }}), +where {\displaystyle \rho ({\boldsymbol {\Sigma }}_{\epsilon })} \rho ({\boldsymbol \Sigma }_{{\epsilon }}) is an inverse-Wishart distribution and {\displaystyle \rho (\mathbf {B} |{\boldsymbol {\Sigma }}_{\epsilon })} \rho ({\mathbf {B}}|{\boldsymbol \Sigma }_{{\epsilon }}) is some form of normal distribution in the matrix {\displaystyle \mathbf {B} } \mathbf {B} . This is accomplished using the vectorization transformation, which converts the likelihood from a function of the matrices {\displaystyle \mathbf {B} ,{\hat {\mathbf {B} }}} {\mathbf {B}},{\hat {{\mathbf {B}}}} to a function of the vectors {\displaystyle {\boldsymbol {\beta }}={\rm {vec}}(\mathbf {B} ),{\hat {\boldsymbol {\beta }}}={\rm {vec}}({\hat {\mathbf {B} }})} {\boldsymbol \beta }={{\rm {vec}}}({\mathbf {B}}),{\hat {{\boldsymbol \beta }}}={{\rm {vec}}}({\hat {{\mathbf {B}}}}). + +Write + +{\displaystyle {\rm {tr}}((\mathbf {B} -{\hat {\mathbf {B} }})^{\rm {T}}\mathbf {X} ^{\rm {T}}\mathbf {X} (\mathbf {B} -{\hat {\mathbf {B} }}){\boldsymbol {\Sigma }}_{\epsilon }^{-1})={\rm {vec}}(\mathbf {B} -{\hat {\mathbf {B} }})^{\rm {T}}{\rm {vec}}(\mathbf {X} ^{\rm {T}}\mathbf {X} (\mathbf {B} -{\hat {\mathbf {B} }}){\boldsymbol {\Sigma }}_{\epsilon }^{-1})} {{\rm {tr}}}(({\mathbf {B}}-{\hat {{\mathbf {B}}}})^{{{\rm {T}}}}{\mathbf {X}}^{{{\rm {T}}}}{\mathbf {X}}({\mathbf {B}}-{\hat {{\mathbf {B}}}}){\boldsymbol \Sigma }_{{\epsilon }}^{{-1}})={{\rm {vec}}}({\mathbf {B}}-{\hat {{\mathbf {B}}}})^{{{\rm {T}}}}{{\rm {vec}}}({\mathbf {X}}^{{{\rm {T}}}}{\mathbf {X}}({\mathbf {B}}-{\hat {{\mathbf {B}}}}){\boldsymbol \Sigma }_{{\epsilon }}^{{-1}}) +Let + +{\displaystyle {\rm {vec}}(\mathbf {X} ^{\rm {T}}\mathbf {X} (\mathbf {B} -{\hat {\mathbf {B} }}){\boldsymbol {\Sigma }}_{\epsilon }^{-1})=({\boldsymbol {\Sigma }}_{\epsilon }^{-1}\otimes \mathbf {X} ^{\rm {T}}\mathbf {X} ){\rm {vec}}(\mathbf {B} -{\hat {\mathbf {B} }}),} {{\rm {vec}}}({\mathbf {X}}^{{{\rm {T}}}}{\mathbf {X}}({\mathbf {B}}-{\hat {{\mathbf {B}}}}){\boldsymbol \Sigma }_{{\epsilon }}^{{-1}})=({\boldsymbol \Sigma }_{{\epsilon }}^{{-1}}\otimes {\mathbf {X}}^{{{\rm {T}}}}{\mathbf {X}}){{\rm {vec}}}({\mathbf {B}}-{\hat {{\mathbf {B}}}}), +where {\displaystyle \mathbf {A} \otimes \mathbf {B} } {\mathbf {A}}\otimes {\mathbf {B}} denotes the Kronecker product of matrices A and B, a generalization of the outer product which multiplies an {\displaystyle m\times n} m\times n matrix by a {\displaystyle p\times q} p\times q matrix to generate an {\displaystyle mp\times nq} mp\times nq matrix, consisting of every combination of products of elements from the two matrices. + +Then + +{\displaystyle {\rm {vec}}(\mathbf {B} -{\hat {\mathbf {B} }})^{\rm {T}}({\boldsymbol {\Sigma }}_{\epsilon }^{-1}\otimes \mathbf {X} ^{\rm {T}}\mathbf {X} ){\rm {vec}}(\mathbf {B} -{\hat {\mathbf {B} }})} {{\rm {vec}}}({\mathbf {B}}-{\hat {{\mathbf {B}}}})^{{{\rm {T}}}}({\boldsymbol \Sigma }_{{\epsilon }}^{{-1}}\otimes {\mathbf {X}}^{{{\rm {T}}}}{\mathbf {X}}){{\rm {vec}}}({\mathbf {B}}-{\hat {{\mathbf {B}}}}) +{\displaystyle =({\boldsymbol {\beta }}-{\hat {\boldsymbol {\beta }}})^{\rm {T}}({\boldsymbol {\Sigma }}_{\epsilon }^{-1}\otimes \mathbf {X} ^{\rm {T}}\mathbf {X} )({\boldsymbol {\beta }}-{\hat {\boldsymbol {\beta }}})} =({\boldsymbol \beta }-{\hat {{\boldsymbol \beta }}})^{{{\rm {T}}}}({\boldsymbol \Sigma }_{{\epsilon }}^{{-1}}\otimes {\mathbf {X}}^{{{\rm {T}}}}{\mathbf {X}})({\boldsymbol \beta }-{\hat {{\boldsymbol \beta }}}) +which will lead to a likelihood which is normal in {\displaystyle ({\boldsymbol {\beta }}-{\hat {\boldsymbol {\beta }}})} ({\boldsymbol \beta }-{\hat {{\boldsymbol \beta }}}). + +With the likelihood in a more tractable form, we can now find a natural (conditional) conjugate prior. + +Conjugate prior distribution[edit] +The natural conjugate prior using the vectorized variable {\displaystyle {\boldsymbol {\beta }}} {\boldsymbol {\beta }} is of the form[1]: + +{\displaystyle \rho ({\boldsymbol {\beta }},{\boldsymbol {\Sigma }}_{\epsilon })=\rho ({\boldsymbol {\Sigma }}_{\epsilon })\rho ({\boldsymbol {\beta }}|{\boldsymbol {\Sigma }}_{\epsilon })} {\displaystyle \rho ({\boldsymbol {\beta }},{\boldsymbol {\Sigma }}_{\epsilon })=\rho ({\boldsymbol {\Sigma }}_{\epsilon })\rho ({\boldsymbol {\beta }}|{\boldsymbol {\Sigma }}_{\epsilon })}, +where + +{\displaystyle \rho ({\boldsymbol {\Sigma }}_{\epsilon })\sim {\mathcal {W}}^{-1}(\mathbf {V_{0}} ,{\boldsymbol {\nu }}_{0})} {\displaystyle \rho ({\boldsymbol {\Sigma }}_{\epsilon })\sim {\mathcal {W}}^{-1}(\mathbf {V_{0}} ,{\boldsymbol {\nu }}_{0})} +and + +{\displaystyle \rho ({\boldsymbol {\beta }}|{\boldsymbol {\Sigma }}_{\epsilon })\sim N({\boldsymbol {\beta }}_{0},{\boldsymbol {\Sigma }}_{\epsilon }\otimes {\boldsymbol {\Lambda }}_{0}).} {\displaystyle \rho ({\boldsymbol {\beta }}|{\boldsymbol {\Sigma }}_{\epsilon })\sim N({\boldsymbol {\beta }}_{0},{\boldsymbol {\Sigma }}_{\epsilon }\otimes {\boldsymbol {\Lambda }}_{0}).} +Posterior distribution[edit] +Using the above prior and likelihood, the posterior distribution can be expressed as[1]: + +{\displaystyle \rho ({\boldsymbol {\beta }},{\boldsymbol {\Sigma }}_{\epsilon }|\mathbf {Y} ,\mathbf {X} )\propto |{\boldsymbol {\Sigma }}_{\epsilon }|^{-({\boldsymbol {\nu }}_{0}+m+1)/2}\exp {(-{\frac {1}{2}}{\rm {tr}}(\mathbf {V_{0}} {\boldsymbol {\Sigma }}_{\epsilon }^{-1}))}} {\displaystyle \rho ({\boldsymbol {\beta }},{\boldsymbol {\Sigma }}_{\epsilon }|\mathbf {Y} ,\mathbf {X} )\propto |{\boldsymbol {\Sigma }}_{\epsilon }|^{-({\boldsymbol {\nu }}_{0}+m+1)/2}\exp {(-{\frac {1}{2}}{\rm {tr}}(\mathbf {V_{0}} {\boldsymbol {\Sigma }}_{\epsilon }^{-1}))}} +{\displaystyle \times |{\boldsymbol {\Sigma }}_{\epsilon }|^{-k/2}\exp {(-{\frac {1}{2}}{\rm {tr}}((\mathbf {B} -\mathbf {B_{0}} )^{\rm {T}}{\boldsymbol {\Lambda }}_{0}(\mathbf {B} -\mathbf {B_{0}} ){\boldsymbol {\Sigma }}_{\epsilon }^{-1}))}} {\displaystyle \times |{\boldsymbol {\Sigma }}_{\epsilon }|^{-k/2}\exp {(-{\frac {1}{2}}{\rm {tr}}((\mathbf {B} -\mathbf {B_{0}} )^{\rm {T}}{\boldsymbol {\Lambda }}_{0}(\mathbf {B} -\mathbf {B_{0}} ){\boldsymbol {\Sigma }}_{\epsilon }^{-1}))}} +{\displaystyle \times |{\boldsymbol {\Sigma }}_{\epsilon }|^{-n/2}\exp {(-{\frac {1}{2}}{\rm {tr}}((\mathbf {Y} -\mathbf {XB} )^{\rm {T}}(\mathbf {Y} -\mathbf {XB} ){\boldsymbol {\Sigma }}_{\epsilon }^{-1}))},} {\displaystyle \times |{\boldsymbol {\Sigma }}_{\epsilon }|^{-n/2}\exp {(-{\frac {1}{2}}{\rm {tr}}((\mathbf {Y} -\mathbf {XB} )^{\rm {T}}(\mathbf {Y} -\mathbf {XB} ){\boldsymbol {\Sigma }}_{\epsilon }^{-1}))},} +where {\displaystyle {\rm {vec}}(\mathbf {B_{0}} )={\boldsymbol {\beta }}_{0}} {\displaystyle {\rm {vec}}(\mathbf {B_{0}} )={\boldsymbol {\beta }}_{0}}. The terms involving {\displaystyle \mathbf {B} } \mathbf {B} can be grouped (with {\displaystyle {\boldsymbol {\Lambda }}_{0}=\mathbf {U} ^{\rm {T}}\mathbf {U} } {\displaystyle {\boldsymbol {\Lambda }}_{0}=\mathbf {U} ^{\rm {T}}\mathbf {U} }) using: + +{\displaystyle (\mathbf {B} -\mathbf {B_{0}} )^{\rm {T}}{\boldsymbol {\Lambda }}_{0}(\mathbf {B} -\mathbf {B_{0}} )+(\mathbf {Y} -\mathbf {XB} )^{\rm {T}}(\mathbf {Y} -\mathbf {XB} )} {\displaystyle (\mathbf {B} -\mathbf {B_{0}} )^{\rm {T}}{\boldsymbol {\Lambda }}_{0}(\mathbf {B} -\mathbf {B_{0}} )+(\mathbf {Y} -\mathbf {XB} )^{\rm {T}}(\mathbf {Y} -\mathbf {XB} )} +{\displaystyle =\left({\begin{bmatrix}\mathbf {Y} \\\mathbf {UB_{0}} \end{bmatrix}}-{\begin{bmatrix}\mathbf {X} \\\mathbf {U} \end{bmatrix}}\mathbf {B} \right)^{\rm {T}}\left({\begin{bmatrix}\mathbf {Y} \\\mathbf {UB_{0}} \end{bmatrix}}-{\begin{bmatrix}\mathbf {X} \\\mathbf {U} \end{bmatrix}}\mathbf {B} \right)} {\displaystyle =\left({\begin{bmatrix}\mathbf {Y} \\\mathbf {UB_{0}} \end{bmatrix}}-{\begin{bmatrix}\mathbf {X} \\\mathbf {U} \end{bmatrix}}\mathbf {B} \right)^{\rm {T}}\left({\begin{bmatrix}\mathbf {Y} \\\mathbf {UB_{0}} \end{bmatrix}}-{\begin{bmatrix}\mathbf {X} \\\mathbf {U} \end{bmatrix}}\mathbf {B} \right)} +{\displaystyle =\left({\begin{bmatrix}\mathbf {Y} \\\mathbf {UB_{0}} \end{bmatrix}}-{\begin{bmatrix}\mathbf {X} \\\mathbf {U} \end{bmatrix}}\mathbf {B_{n}} \right)^{\rm {T}}\left({\begin{bmatrix}\mathbf {Y} \\\mathbf {UB_{0}} \end{bmatrix}}-{\begin{bmatrix}\mathbf {X} \\\mathbf {U} \end{bmatrix}}\mathbf {B_{n}} \right)+(\mathbf {B} -\mathbf {B_{n}} )^{\rm {T}}(\mathbf {X} ^{\rm {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0})(\mathbf {B} -\mathbf {B_{n}} )} {\displaystyle =\left({\begin{bmatrix}\mathbf {Y} \\\mathbf {UB_{0}} \end{bmatrix}}-{\begin{bmatrix}\mathbf {X} \\\mathbf {U} \end{bmatrix}}\mathbf {B_{n}} \right)^{\rm {T}}\left({\begin{bmatrix}\mathbf {Y} \\\mathbf {UB_{0}} \end{bmatrix}}-{\begin{bmatrix}\mathbf {X} \\\mathbf {U} \end{bmatrix}}\mathbf {B_{n}} \right)+(\mathbf {B} -\mathbf {B_{n}} )^{\rm {T}}(\mathbf {X} ^{\rm {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0})(\mathbf {B} -\mathbf {B_{n}} )} +{\displaystyle =(\mathbf {Y} -\mathbf {XB_{n}} )^{\rm {T}}(\mathbf {Y} -\mathbf {XB_{n}} )+(\mathbf {B_{0}} -\mathbf {B_{n}} )^{\rm {T}}{\boldsymbol {\Lambda }}_{0}(\mathbf {B_{0}} -\mathbf {B_{n}} )+(\mathbf {B} -\mathbf {B_{n}} )^{\rm {T}}(\mathbf {X} ^{\rm {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0})(\mathbf {B} -\mathbf {B_{n}} )} {\displaystyle =(\mathbf {Y} -\mathbf {XB_{n}} )^{\rm {T}}(\mathbf {Y} -\mathbf {XB_{n}} )+(\mathbf {B_{0}} -\mathbf {B_{n}} )^{\rm {T}}{\boldsymbol {\Lambda }}_{0}(\mathbf {B_{0}} -\mathbf {B_{n}} )+(\mathbf {B} -\mathbf {B_{n}} )^{\rm {T}}(\mathbf {X} ^{\rm {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0})(\mathbf {B} -\mathbf {B_{n}} )}, +with + +{\displaystyle \mathbf {B_{n}} =(\mathbf {X} ^{\rm {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0})^{-1}(\mathbf {X} ^{\rm {T}}\mathbf {X} {\hat {\mathbf {B} }}+{\boldsymbol {\Lambda }}_{0}\mathbf {B_{0}} )=(\mathbf {X} ^{\rm {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0})^{-1}(\mathbf {X} ^{\rm {T}}\mathbf {Y} +{\boldsymbol {\Lambda }}_{0}\mathbf {B_{0}} )} {\displaystyle \mathbf {B_{n}} =(\mathbf {X} ^{\rm {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0})^{-1}(\mathbf {X} ^{\rm {T}}\mathbf {X} {\hat {\mathbf {B} }}+{\boldsymbol {\Lambda }}_{0}\mathbf {B_{0}} )=(\mathbf {X} ^{\rm {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0})^{-1}(\mathbf {X} ^{\rm {T}}\mathbf {Y} +{\boldsymbol {\Lambda }}_{0}\mathbf {B_{0}} )}. +This now allows us to write the posterior in a more useful form: + +{\displaystyle \rho ({\boldsymbol {\beta }},{\boldsymbol {\Sigma }}_{\epsilon }|\mathbf {Y} ,\mathbf {X} )\propto |{\boldsymbol {\Sigma }}_{\epsilon }|^{-({\boldsymbol {\nu }}_{0}+m+n+1)/2}\exp {(-{\frac {1}{2}}{\rm {tr}}((\mathbf {V_{0}} +(\mathbf {Y} -\mathbf {XB_{n}} )^{\rm {T}}(\mathbf {Y} -\mathbf {XB_{n}} )+(\mathbf {B_{n}} -\mathbf {B_{0}} )^{\rm {T}}{\boldsymbol {\Lambda }}_{0}(\mathbf {B_{n}} -\mathbf {B_{0}} )){\boldsymbol {\Sigma }}_{\epsilon }^{-1}))}} {\displaystyle \rho ({\boldsymbol {\beta }},{\boldsymbol {\Sigma }}_{\epsilon }|\mathbf {Y} ,\mathbf {X} )\propto |{\boldsymbol {\Sigma }}_{\epsilon }|^{-({\boldsymbol {\nu }}_{0}+m+n+1)/2}\exp {(-{\frac {1}{2}}{\rm {tr}}((\mathbf {V_{0}} +(\mathbf {Y} -\mathbf {XB_{n}} )^{\rm {T}}(\mathbf {Y} -\mathbf {XB_{n}} )+(\mathbf {B_{n}} -\mathbf {B_{0}} )^{\rm {T}}{\boldsymbol {\Lambda }}_{0}(\mathbf {B_{n}} -\mathbf {B_{0}} )){\boldsymbol {\Sigma }}_{\epsilon }^{-1}))}} +{\displaystyle \times |{\boldsymbol {\Sigma }}_{\epsilon }|^{-k/2}\exp {(-{\frac {1}{2}}{\rm {tr}}((\mathbf {B} -\mathbf {B_{n}} )^{\rm {T}}{\boldsymbol {\Lambda }}_{0})(\mathbf {B} -\mathbf {B_{n}} ){\boldsymbol {\Sigma }}_{\epsilon }^{-1}))}} {\displaystyle \times |{\boldsymbol {\Sigma }}_{\epsilon }|^{-k/2}\exp {(-{\frac {1}{2}}{\rm {tr}}((\mathbf {B} -\mathbf {B_{n}} )^{\rm {T}}{\boldsymbol {\Lambda }}_{0})(\mathbf {B} -\mathbf {B_{n}} ){\boldsymbol {\Sigma }}_{\epsilon }^{-1}))}}. +This takes the form of an inverse-Wishart distribution times a Matrix normal distribution: + +{\displaystyle \rho ({\boldsymbol {\Sigma }}_{\epsilon }|\mathbf {Y} ,\mathbf {X} )\sim {\mathcal {W}}^{-1}(\mathbf {V_{n}} ,{\boldsymbol {\nu }}_{n})} {\displaystyle \rho ({\boldsymbol {\Sigma }}_{\epsilon }|\mathbf {Y} ,\mathbf {X} )\sim {\mathcal {W}}^{-1}(\mathbf {V_{n}} ,{\boldsymbol {\nu }}_{n})} +and + +{\displaystyle \rho ({\boldsymbol {\beta }}|\mathbf {Y} ,\mathbf {X} ,{\boldsymbol {\Sigma }}_{\epsilon })\sim {\mathcal {MN}}_{k,m}(\mathbf {B_{n}} ,{\boldsymbol {\Lambda }}_{n}^{-1},{\boldsymbol {\Sigma }}_{\epsilon })} {\displaystyle \rho ({\boldsymbol {\beta }}|\mathbf {Y} ,\mathbf {X} ,{\boldsymbol {\Sigma }}_{\epsilon })\sim {\mathcal {MN}}_{k,m}(\mathbf {B_{n}} ,{\boldsymbol {\Lambda }}_{n}^{-1},{\boldsymbol {\Sigma }}_{\epsilon })}. +The parameters of this posterior are given by: + +{\displaystyle \mathbf {V_{n}} =\mathbf {V_{0}} +(\mathbf {Y} -\mathbf {XB_{n}} )^{\rm {T}}(\mathbf {Y} -\mathbf {XB_{n}} )+(\mathbf {B_{n}} -\mathbf {B_{0}} )^{\rm {T}}{\boldsymbol {\Lambda }}_{0}(\mathbf {B_{n}} -\mathbf {B_{0}} )} {\displaystyle \mathbf {V_{n}} =\mathbf {V_{0}} +(\mathbf {Y} -\mathbf {XB_{n}} )^{\rm {T}}(\mathbf {Y} -\mathbf {XB_{n}} )+(\mathbf {B_{n}} -\mathbf {B_{0}} )^{\rm {T}}{\boldsymbol {\Lambda }}_{0}(\mathbf {B_{n}} -\mathbf {B_{0}} )} +{\displaystyle {\boldsymbol {\nu }}_{n}={\boldsymbol {\nu }}_{0}+n} {\displaystyle {\boldsymbol {\nu }}_{n}={\boldsymbol {\nu }}_{0}+n} +{\displaystyle \mathbf {B_{n}} =(\mathbf {X} ^{\rm {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0})^{-1}(\mathbf {X} ^{\rm {T}}\mathbf {Y} +{\boldsymbol {\Lambda }}_{0}\mathbf {B_{0}} )} {\displaystyle \mathbf {B_{n}} =(\mathbf {X} ^{\rm {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0})^{-1}(\mathbf {X} ^{\rm {T}}\mathbf {Y} +{\boldsymbol {\Lambda }}_{0}\mathbf {B_{0}} )} +{\displaystyle {\boldsymbol {\Lambda }}_{n}=\mathbf {X} ^{\rm {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0}} {\displaystyle {\boldsymbol {\Lambda }}_{n}=\mathbf {X} ^{\rm {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0}} +~