Chapter 1 Basics of Bayesian linear regression

1.1 Bayes’ theorem

Theorem 1.1 (Bayes' theorem) For events $U, K$ and $P(K) \neq 0$, we have

$$P(U \mid K) = \frac{P(K \mid U) P(U)}{P(K)}$$

We denote $U$ as unknown parameters and $K$ as known parameters. We call $P(U)$ prior and $P(K|U)$ likelihood. The Bayes’ theorem gives us the posterior distribution of unknown parameters given the known parameters $$P(U \mid K) \propto P(U) \cdot P(K \mid U)$$

1.2 Normal-Inverse-Gamma (NIG) prior

1.2.1 Joint distribution of NIG prior

Definition 1.1 (Normal-Inverse-Gamma Distribution) Suppose $$\begin{aligned} & \beta \mid \sigma^2, \mu, M \sim N(\mu,\sigma^2 M) \\ & \sigma^2 \mid a, b \sim IG(a, b) \end{aligned}$$ Then $(\beta,\sigma^2)$ has a Normal-Inverse-Gamma distribution, denoted as $(\beta,\sigma^2) \sim NIG(\mu,M,a,b)$.

We use a Normal-Inverse-Gamma prior for $(\beta, \sigma^2)$

$$\begin{align} P(\beta, \sigma^{2}) &= NIG \left(\beta, \sigma^{2} \mid m_{0}, M_{0}, a_{0}, b_{0}\right) \\ &= \frac{b_0^{a_0}}{\Gamma\left(a_{0}\right)} \left(\frac{1}{\sigma^{2}}\right)^{a_{0}+1} e^{-\frac{b_{0}}{\sigma^{2}}} \frac{1}{(2 \pi \sigma^{2})^{\frac{p}{2}}\left| M_{0}\right|^{\frac{1}{2}}} e^{-\frac{1}{2 \sigma^{2}} Q \left(\beta, m_{0}, M_{0}\right)} \end{align}$$

where $Q(x, m, M)=(x-m)^{\mathrm{\scriptscriptstyle T}} M^{-1} (x-m) \, .$

$$\begin{align} P(\beta, \sigma^{2}) &= NIG \left(\beta, \sigma^{2} \mid m_{0}, M_{0}, a_{0}, b_{0}\right) \\ &= IG\left(\sigma^{2} \mid a_{0}, b_{0}\right) \cdot N\left(\beta \mid m_{0}, \sigma^{2} M_{0}\right) \\ &= \frac{b_0^{a_0}}{\Gamma\left(a_{0}\right)} \left(\frac{1}{\sigma^{2}}\right)^{a_{0}+1} e^{-\frac{b_{0}}{\sigma^{2}}} \frac{1}{(2 \pi)^{\frac{p}{2}}\left|\sigma^{2} M_{0}\right|^{\frac{1}{2}}} e^{-\frac{1}{2 \sigma^{2}} Q \left(\beta, m_{0}, M_{0}\right)} \\ &= \frac{b_0^{a_0}}{\Gamma\left(a_{0}\right)} \left(\frac{1}{\sigma^{2}}\right)^{a_{0}+1} e^{-\frac{b_{0}}{\sigma^{2}}} \frac{1}{(2 \pi \sigma^{2})^{\frac{p}{2}}\left| M_{0}\right|^{\frac{1}{2}}} e^{-\frac{1}{2 \sigma^{2}} Q \left(\beta, m_{0}, M_{0}\right)} \\ \end{align}$$

where $Q(x, m, M)=(x-m)^{\mathrm{\scriptscriptstyle T}} M^{-1} (x-m) \, .$

Note: the Inverse-Gamma ($IG$) distribution has a relationship with Gamma distribution. Given $X \sim Gamma(\alpha, \beta)$, the density function of $X$ is $f(x)=\frac{\beta^{\alpha}}{\Gamma(\alpha)} x^{\alpha-1} e^{-\beta x}$. Then $Y=\frac{1}{X} \sim IG\left(\alpha, \beta\right)$ with density function $f(y)=\frac{\beta^{\alpha}}{\Gamma(\alpha)} x^{-\alpha-1} e^{-\frac{\beta}{x}}$.

1.2.2 Marginal distribution of NIG prior

As for marginal priors, we can can get it by integration

$$\begin{aligned} P(\sigma^2) & = \int N I G\left(\beta, \sigma^{2} \mid m_{0}, M_{0}, a_{0}, b_{0}\right) \ d\beta=I G\left(\sigma^{2} \mid a_{0}, b_{0}\right) \\ P(\beta) & = \int N I G\left(\beta, \sigma^{2} \mid m_{0}, M_{0}, a_{0}, b_{0}\right) \ d\sigma^{2}=t_{2a_0}\left(m_0, \frac{b_0}{a_0}M_0\right) \\ \end{aligned}$$

$$\begin{align} P\left(\sigma^{2} \right) &= \int NIG\left(\beta, \sigma^{2} \mid m_{0}, M_{0}, a_{0}, b_{0}\right) \ d\beta \\ &=IG\left(\sigma^{2} \mid a_{0}, b_{0}\right) \int N\left(\beta \mid m_{0}, \sigma^{2} M_{0}\right) \ d\beta \\ &=IG\left(\sigma^{2} \mid a_{0}, b_{0}\right) \end{align}$$

$$\begin{align} P(\beta ) &=\int NIG \left(\beta, \sigma^{2} \mid m_{0}, M_{0}, a_{0}, b_{0}\right) \ d\sigma^{2} \\ &= \int \frac{b_0^{a_0}}{\Gamma\left(a_{0}\right)} \left(\frac{1}{\sigma^{2}}\right)^{a_{0}+1} e^{-\frac{b_{0}}{\sigma^{2}}} \frac{1}{(2 \pi \sigma^{2})^{\frac{p}{2}}\left| M_{0}\right|^{\frac{1}{2}}} e^{-\frac{1}{2 \sigma^{2}} Q \left(\beta, m_{0}, M_{0}\right)} \ d\sigma^{2} \\ &= \frac{b_0^{a_0}}{\Gamma\left(a_{0}\right)(2 \pi)^{\frac{p}{2}}\left| M_{0}\right|^{\frac{1}{2}}} \int \left(\frac{1}{\sigma^{2}}\right)^{a_{0}+\frac{p}{2}+1} e^{-\frac{1}{\sigma^{2}}\left(b_{0}+\frac{1}{2} Q \left(\beta, m_{0}, M_{0}\right)\right)} \ d\sigma^{2} \\ & \quad \left(\text{let } u = \frac{1}{\sigma^2}, \left|d\sigma^{2}\right|=\frac{1}{u^2} d u\right) \\ &= \frac{b_0^{a_0}}{\Gamma\left(a_{0}\right)(2 \pi)^{\frac{p}{2}}\left| M_{0}\right|^{\frac{1}{2}}} \int u^{a_{0}+\frac{p}{2}+1} e^{-\left(b_{0}+\frac{1}{2} Q \left(\beta, m_{0}, M_{0}\right)\right) u } \frac{1}{u^2} \ du \\ &= \frac{b_0^{a_0}}{\Gamma\left(a_{0}\right)(2 \pi)^{\frac{p}{2}}\left| M_{0}\right|^{\frac{1}{2}}} \int u^{a_{0}+\frac{p}{2}-1} e^{-\left(b_{0}+\frac{1}{2} Q \left(\beta, m_{0}, M_{0}\right)\right) u} \ du \\ & \quad \left(\text{by Gamma integral function:} \int x^{\alpha - 1} exp^{-\beta x} dx = \frac{\Gamma(\alpha)}{\beta^{\alpha}}\right) \\ &= \frac{b_{0}^{a_{0}} }{\Gamma\left(a_{0}\right)(2 \pi)^\frac{p}{2}\left|M_{0}\right|^{\frac{1}{2}}} \frac{\Gamma\left(a_{0}+\frac{p}{2}\right)}{\left(b_{0}+\frac{1}{2} Q(\beta,m_0,M_0)\right)^{\left(a_{0}+\frac{p}{2}\right)}} \\ & = \frac{b_0^{a_0}\Gamma\left(a_{0}+\frac{p}{2}\right)}{\Gamma\left(a_{0}\right)(2 \pi)^ \frac{p}{2}\left|M_{0}\right|^{\frac{1}{2}}} \left(b_0(1+\frac{1}{2 b_0} Q(\beta,m_0,M_0))\right)^{-\left(a_{0}+\frac{p}{2}\right)} \\ & = \frac{b_0^{a_0}\Gamma\left(a_{0}+\frac{p}{2}\right) b_0^{- \left( a_0+\frac{p}{2}\right)}}{\Gamma\left(a_{0}\right)(2 \pi)^ \frac{p}{2}\left|M_{0}\right|^{\frac{1}{2}}} \left(1+\frac{1}{2 b_0} \left(\beta-m_{0}\right)^{\mathrm{\scriptscriptstyle T}} M_{0}^{-1}\left(\beta-m_{0}\right) \right)^{-\left(a_{0}+\frac{p}{2}\right)} \\ & =\frac{\Gamma\left(a_{0}+\frac{p}{2}\right)}{\left(2 \pi \right)^{\frac{p}{2}} b_{0}^{\frac{p}{2}} \Gamma\left(a_{0}\right)|M|^{\frac{1}{2}}}\left(1+\frac{1}{2 b_{0}} \left(\beta-m_{0}\right)^{\mathrm{\scriptscriptstyle T}} M_{0}^{-1}\left(\beta-m_{0}\right) \right)^{-\left(a_{0}+\frac{p}{2}\right)} \\ & =\frac{\Gamma\left(a_{0}+\frac{p}{2}\right)} {\left(2 \pi \right)^{\frac{p}{2}}\left(a_{0} \cdot \frac{b_{0}}{a_{0}}\right)^{\frac{p}{2}} \Gamma\left(a_{0}\right)|M|^{\frac{1}{2}}} \left(1+\frac{1}{2 a_{0} \cdot \frac{b_{0}}{a_{0}}} \left(\beta-m_{0}\right)^{\mathrm{\scriptscriptstyle T}} M_{0}^{-1}\left(\beta-m_{0}\right)\right)^{-\left(a_{0}+\frac{p}{2}\right)}\\ & =\frac{\Gamma\left(a_{0}+\frac{p}{2}\right)}{\left(2 a_{0} \pi\right)^{\frac{p}{2}} \Gamma\left(a_{0}\right)\left|\frac{b_{0}}{a_{0}} M\right|^{\frac{1}{2}}}\left(1+\frac{1}{2 a_{0}} \left(\beta-m_{0}\right)^{\mathrm{\scriptscriptstyle T}}\left(\frac{b_{0}}{a_{0}} M_{0}\right)^{-1}\left(\beta-m_{0}\right)\right)^{-\left(a_{0}+\frac{p}{2}\right)} \\ & =t_{2a_0}\left(m_0, \frac{b_0}{a_0}M_0\right) \; \end{align}$$

Note: the density of multivariate t-distribution is given by

$$t_v(\mu, \Sigma)=\frac{\Gamma\left(\frac{v+p}{2}\right)}{(v \pi)^{\frac{p}{2}} \Gamma\left(\frac{v}{2}\right) |\Sigma|^{\frac{1}{2}}}\left(1+\frac{1}{v}(x-\mu)^{\mathrm{\scriptscriptstyle T}} \Sigma^{-1}(x-\mu)\right)^{-\frac{v+p}{2}}$$

1.3 Conjugate Bayesian linear regression and M&m formula

Let $y_{n \times 1}$ be outcome variable and $X_{n \times p}$ be corresponding covariates. Assume $V$ is known. The model is given by

$$\begin{align} & y=X \beta+\epsilon \ , \ \epsilon \sim N\left(0, \sigma^2 V\right) \\ & \beta=m_0+\omega \ , \ \omega \sim N\left(0, \sigma^2 M_0\right) \\ & \sigma^2 \sim I G\left(a_0, b_0\right) \end{align}$$

The posterior distribution of $(\beta, \sigma^2)$ is given by

$$\begin{align} P\left(\beta, \sigma^{2} \mid y\right) = NIG\left(\beta, \sigma^{2} \mid M_{1}m_{1}, M_{1}, a_{1}, b_{1}\right) \; \end{align}$$

where

$$\begin{align} M_{1}^{-1} &= M_{0}^{-1}+X^{\mathrm{\scriptscriptstyle T}} V^{-1} X \\ m_{1}&=M_{0}^{-1} m_{0}+X^{\mathrm{\scriptscriptstyle T}} V^{-1} y \\ a_{1}&=a_{0}+\frac{p}{2} \\ b_{1}&=b_{0}+\frac{c^{\ast}}{2}= b_{0}+\frac{1}{2}\left(m_{0}^{\mathrm{\scriptscriptstyle T}} M_{0}^{-1} m_{0}+y^{\mathrm{\scriptscriptstyle T}} V^{-1} y-m_{1}^{\mathrm{\scriptscriptstyle T}} M_{1} m_{1}\right) \end{align}$$

$$\begin{align}\label{eq:post_dist} P\left(\beta, \sigma^{2} \mid y\right) & \propto NIG\left(\beta, \sigma^{2} \mid m_{0}, M_{0}, a_{0}, b_{0}\right) \cdot N\left(y \mid X \beta, \sigma^{2} V\right) \nonumber\\ & \propto IG\left(\sigma^{2} \mid a_{0}, b_{0}\right) \cdot N\left(\beta \mid m_{0}, \sigma^{2} M_{0}\right) \cdot N\left(y \mid X \beta, \sigma^{2} V\right) \nonumber\\ & \propto \frac{b_0^{a_0}}{\Gamma\left(a_{0}\right)} \left(\frac{1}{\sigma^{2}}\right)^{a_{0}+1} e^{-\frac{b_{0}}{\sigma^{2}}} \frac{1}{(2 \pi \sigma^{2})^{\frac{p}{2}}\left| M_{0}\right|^{\frac{1}{2}}} e^{-\frac{1}{2 \sigma^{2}} Q \left(\beta, m_{0}, M_{0}\right)} \frac{1}{(2 \pi \sigma^{2})^{\frac{p}{2}}\left| V\right|^{\frac{1}{2}}} e^{-\frac{1}{2 \sigma^{2}} Q \left(y, X \beta, V\right)} \nonumber\\ & \propto \left(\frac{1}{\sigma^{2}}\right)^{a_{0}+p+1} e^{-\frac{b_{0}}{\sigma^{2}}} e^{-\frac{1}{2 \sigma^{2}} \left(Q \left(\beta, m_{0}, M_{0}\right)+Q \left(y, X \beta, V\right)\right)}\; \end{align}$$

where

$$\begin{align}\label{eq:multivariate_completion_square} Q \left(\beta, m_{0}, M_{0}\right)+Q \left(y, X \beta, V\right) &= (\beta - m_{0})^{\mathrm{\scriptscriptstyle T}}M_{0}^{-1}(\beta - m_{0}) + (y - X\beta)^{\mathrm{\scriptscriptstyle T}}V^{-1}(y - X\beta)\; \nonumber\\ &= \beta^{\mathrm{\scriptscriptstyle T}}M_{0}^{-1}\beta - 2\beta^{\mathrm{\scriptscriptstyle T}}M_{0}^{-1}m_{0} + m_{0}^{\mathrm{\scriptscriptstyle T}}M_{0}^{-1}m_{0} \nonumber\\ &\qquad + \beta^{\mathrm{\scriptscriptstyle T}}X^{\mathrm{\scriptscriptstyle T}}V^{-1}X\beta - 2\beta^{\mathrm{\scriptscriptstyle T}} X^{\mathrm{\scriptscriptstyle T}}V^{-1}y + y^{\mathrm{\scriptscriptstyle T}}V^{-1}y \nonumber\\ &= \beta^{\mathrm{\scriptscriptstyle T}} \left(M_{0}^{-1} + X^{\mathrm{\scriptscriptstyle T}}V^{-1}X\right) \beta - 2\beta^{\mathrm{\scriptscriptstyle T}}\left(M_{0}^{-1}m_{0} + X^{\mathrm{\scriptscriptstyle T}}V^{-1}y\right) \nonumber\\ &\qquad + m_{0}^{\mathrm{\scriptscriptstyle T}} M_{0}^{-1}m_{0} + y^{\mathrm{\scriptscriptstyle T}}V^{-1}y \nonumber \\ &= \beta^{\mathrm{\scriptscriptstyle T}}M_{1}^{-1}\beta - 2\beta^{\mathrm{\scriptscriptstyle T}} m_{1} + c\nonumber\\ &= (\beta - M_{1}m_{1})^{\mathrm{\scriptscriptstyle T}}M_{1}^{-1}(\beta - M_{1}m_{1}) - m_{1}^{\mathrm{\scriptscriptstyle T}}M_{1}m_{1} +c \nonumber\\ &= (\beta - M_{1}m_{1})^{\mathrm{\scriptscriptstyle T}}M_{1}^{-1}(\beta - M_{1}m_{1}) +c^{\ast}\; \end{align}$$

where $M_{1}$ is a symmetric positive definite matrix, $m_{1}$ is a vector, and $c$ & $c^{\ast}$ are scalars given by

$$\begin{align} M_{1}^{-1} &= M_{0}^{-1} + X^{\mathrm{\scriptscriptstyle T}}V^{-1}X \\ m_{1} &= M_{0}^{-1}m_{0} + X^{\mathrm{\scriptscriptstyle T}}V^{-1}y \\ c &= m_{0}^{\mathrm{\scriptscriptstyle T}} M_{0}^{-1}m_{0} + y^{\mathrm{\scriptscriptstyle T}}V^{-1}y \\ c^{\ast} &= c - m^{\mathrm{\scriptscriptstyle T}}Mm = m_{0}^{\mathrm{\scriptscriptstyle T}} M_{0}^{-1}m_{0} + y^{\mathrm{\scriptscriptstyle T}}V^{-1}y - m_{1}^{\mathrm{\scriptscriptstyle T}}M_{1}m_{1} \end{align}$$

Note: $M_{1}$, $m_{1}$ and $c$ do not depend upon $\beta$.

Then, we have

$$\begin{align} P\left(\beta, \sigma^{2} \mid y\right) & \propto \left(\frac{1}{\sigma^{2}}\right)^{a_{0}+p+1} e^{-\frac{b_{0}}{\sigma^{2}}} e^{-\frac{1}{2 \sigma^{2}} ((\beta - M_{1}m_{1})^{\mathrm{\scriptscriptstyle T}}M_{1}^{-1}(\beta - M_{1}m_{1}) +c^{\ast})}\\ & \propto \left(\frac{1}{\sigma^{2}}\right)^{a_{0}+p+1} e^{-\frac{b_{0}+\frac{c^{\ast}}{2}}{\sigma^{2}}} e^{-\frac{1}{2 \sigma^{2}} (\beta - M_{1}m_{1})^{\mathrm{\scriptscriptstyle T}}M_{1}^{-1}(\beta - M_{1}m_{1})}\\ & \propto \left(\frac{1}{\sigma^{2}}\right)^{a_{0}+\frac{p}{2}+1} e^{-\frac{b_{0}+\frac{c^{\ast}}{2}}{\sigma^{2}}} (\frac{1}{\sigma^2})^{\frac{p}{2}} e^{-\frac{1}{2 \sigma^{2}} (\beta - M_{1}m_{1})^{\mathrm{\scriptscriptstyle T}}M_{1}^{-1}(\beta - M_{1}m_{1})}\\ & \propto IG\left(\sigma^{2} \mid a_{0}+\frac{p}{2}, b_{0}+\frac{c^{\ast}}{2} \right) \cdot N\left(\beta \mid M_{1}m_{1}, \sigma^{2} M_{1}\right) \\ & \propto IG\left(\sigma^{2} \mid a_{1}, b_{1} \right) \cdot N\left(\beta \mid M_{1}m_{1}, \sigma^{2} M_{1}\right) \\ & \propto NIG\left(\beta, \sigma^{2} \mid M_{1}m_{1}, M_{1}, a_{1}, b_{1}\right) \; \end{align}$$

where

$$\begin{align} M_{1}^{-1} &= M_{0}^{-1}+X^{\mathrm{\scriptscriptstyle T}} V^{-1} X \\ m_{1}&=M_{0}^{-1} m_{0}+X^{\mathrm{\scriptscriptstyle T}} V^{-1} y \\ a_{1}&=a_{0}+\frac{p}{2} \\ b_{1}&=b_{0}+\frac{c^{\ast}}{2}= b_{0}+\frac{1}{2}\left(m_{0}^{\mathrm{\scriptscriptstyle T}} M_{0}^{-1} m_{0}+y^{\mathrm{\scriptscriptstyle T}} V^{-1} y-m_{1}^{\mathrm{\scriptscriptstyle T}} M_{1} m_{1}\right) \end{align}$$

From derivation in marginal priors, the marginal posterior distributions can be easily get by updating corresponding parameters

$$\begin{aligned} &P\left(\sigma^{2} \mid y\right)=I G\left(\sigma^{2} \mid a_{1}, b_{1}\right) \\ &P(\beta \mid y)=t_{2a_1}\left(M_1m_1, \frac{b_1}{a_1}M_1\right) \end{aligned}$$

1.4 Updating form of the posterior distribution

We will use two ways to derive the updating form of the posterior distribution.

1.4.1 Method 1: Sherman-Woodbury-Morrison identity

Theorem 1.2 (Sherman-Woodbury-Morrison identity) We have $$\begin{equation}\label{ShermanWoodburyMorrison} \left(A + BDC\right)^{-1} = A^{-1} - A^{-1}B\left(D^{-1}+CA^{-1}B\right)^{-1}CA^{-1} \end{equation}$$ where $A$ and $D$ are square matrices that are invertible and $B$ and $C$ are rectangular (square if $A$ and $D$ have the same dimensions) matrices such that the multiplications are well-defined.

Sherman-Woodbury-Morrison identity is easily verified by multiplying the right hand side with $A + BDC$ and simplifying to reduce it to the identity matrix. Using this formula, we have

$$\begin{aligned} M_1 & = (M_{0}^{-1} + X^{\mathrm{\scriptscriptstyle T}}V^{-1}X)^{-1} \\ & = M_0-M_0 X^{\mathrm{\scriptscriptstyle T}}\left(V+X M_0 X^{\mathrm{\scriptscriptstyle T}}\right)^{-1} X M_0 \\ & = M_0-M_0 X^{\mathrm{\scriptscriptstyle T}} Q^{-1} X M_0 \end{aligned}$$

where $Q = V + X M_0 X^{\mathrm{\scriptscriptstyle T}}$

We can show that $$\begin{align} M_1 m_1 & =m_0+M_0 X^{\mathrm{\scriptscriptstyle T}} Q^{-1}\left(y-X m_0\right) \end{align}$$

$$\begin{align} M_1 m_1 & = \left(M_0^{-1}+X^{\mathrm{\scriptscriptstyle T}} V^{-1} X\right)^{-1} m_1 \\ & = \left(M_0-M_0 X^{\mathrm{\scriptscriptstyle T}}\left(V+X M_0 X^{\mathrm{\scriptscriptstyle T}}\right)^{-1} X M_0\right)m_1 \\ & = \left(M_0-M_0 X^{\mathrm{\scriptscriptstyle T}} Q^{-1} X M_0\right) m_1 \\ & = \left(M_0-M_0 X^{\mathrm{\scriptscriptstyle T}} Q^{-1} X M_0\right)\left(M_0^{-1} m_0+X^{\mathrm{\scriptscriptstyle T}} V^{-1} y\right) \\ & = m_0+M_0 X^{\mathrm{\scriptscriptstyle T}} V^{-1} y-M_0 X^{\mathrm{\scriptscriptstyle T}} Q^{-1} X m_0 - M_0 X^{\mathrm{\scriptscriptstyle T}} Q^{-1} X M_0 X^{\mathrm{\scriptscriptstyle T}} V^{-1} y \\ & = m_0+M_0 X^{\mathrm{\scriptscriptstyle T}}\left(I-Q^{-1} X M_0 X^{\mathrm{\scriptscriptstyle T}}\right) V^{-1} y - M_0 X^{\mathrm{\scriptscriptstyle T}} Q^{-1} X m_0 \\ & = m_0+M_0 X^{\mathrm{\scriptscriptstyle T}} Q^{-1}\left(Q-X M_0 X^{\mathrm{\scriptscriptstyle T}}\right)V^{-1} y - M_0 X^{\mathrm{\scriptscriptstyle T}} Q^{-1} X m_0 \\ & \left(\text { since } Q=V+X M_0 X^{\mathrm{\scriptscriptstyle T}}\right) \\ & = m_0+M_0 X^{\mathrm{\scriptscriptstyle T}} Q^{-1}V V^{-1} y-M_0 X^{\mathrm{\scriptscriptstyle T}} Q^{-1} X m_0 \\ & = m_0+M_0 X^{\mathrm{\scriptscriptstyle T}} Q^{-1} y-M_0 X^{\mathrm{\scriptscriptstyle T}} Q^{-1} X m_0 \\ & = m_0+M_0 X^{\mathrm{\scriptscriptstyle T}} Q^{-1}\left(y-X m_0\right) \\ \end{align}$$

Furthermore, we can simplify that $$\begin{align} m_0^{\mathrm{\scriptscriptstyle T}} M_0^{-1} m_0+y^{\mathrm{\scriptscriptstyle T}} V^{-1} y-m_1^{\mathrm{\scriptscriptstyle T}} M_1 m_1 & = \left(y-X m_0\right)^{\mathrm{\scriptscriptstyle T}} Q^{-1}\left(y-X m_0\right) \end{align}$$

$$\begin{align} & \quad \ m_0^{\mathrm{\scriptscriptstyle T}} M_0^{-1} m_0+y^{\mathrm{\scriptscriptstyle T}} V^{-1} y-m_1^{\mathrm{\scriptscriptstyle T}} M_1 m_1 \\ & = m_0^{\mathrm{\scriptscriptstyle T}} M_0^{-1} m_0+y^{\mathrm{\scriptscriptstyle T}} V^{-1} y-m_1^{\mathrm{\scriptscriptstyle T}} \left(m_0+M_0 X^{\mathrm{\scriptscriptstyle T}} Q^{-1} (y - X m_0)\right) \\ & = m_0^{\mathrm{\scriptscriptstyle T}} M_0^{-1} m_0+y^{\mathrm{\scriptscriptstyle T}} V^{-1} y-m_1^{\mathrm{\scriptscriptstyle T}} m_0 - m_1^{\mathrm{\scriptscriptstyle T}} M_0 X^{\mathrm{\scriptscriptstyle T}} Q^{-1}\left(y-X m_0\right) \\ & = m_0^{\mathrm{\scriptscriptstyle T}} M_0^{-1} m_0+y^{\mathrm{\scriptscriptstyle T}} V^{-1} y -m_0^{\mathrm{\scriptscriptstyle T}}\left(M_0^{-1} m_0+X^{\mathrm{\scriptscriptstyle T}} V^{-1} y\right) \\ & \qquad \qquad \qquad - m_1^{\mathrm{\scriptscriptstyle T}} M_0 X^{\mathrm{\scriptscriptstyle T}} Q^{-1}\left(y-X m_0\right) \\ & = y^{\mathrm{\scriptscriptstyle T}} V^{-1} y-y^{\mathrm{\scriptscriptstyle T}} V^{-1} X m_0 - m_1^{\mathrm{\scriptscriptstyle T}} M_0 X^{\mathrm{\scriptscriptstyle T}} Q^{-1}\left(y-X m_0\right) \\ & = y^{\mathrm{\scriptscriptstyle T}} V^{-1}\left(y-X m_0 \right)-m_1^{\mathrm{\scriptscriptstyle T}} M_0 X^{\mathrm{\scriptscriptstyle T}} Q^{-1}\left(y-X m_0\right) \\ & =y^{\mathrm{\scriptscriptstyle T}} V^{-1}\left(y-X m_0\right)-\underbrace{m_1^{\mathrm{\scriptscriptstyle T}} M_0 X^{\mathrm{\scriptscriptstyle T}} Q^{-1}\left(y-X m_0\right)}_{\substack{\text { simplify from left to right }}} \\ & =y^{\mathrm{\scriptscriptstyle T}} V^{-1}\left(y-X m_0\right)-\left(M_0 m_1\right)^{\mathrm{\scriptscriptstyle T}} X^{\mathrm{\scriptscriptstyle T}} Q^{-1}\left(y-X m_0\right) \\ & =y^{\mathrm{\scriptscriptstyle T}} V^{-1}\left(y-X m_0\right)-\left(m_0+M_0 X^{\mathrm{\scriptscriptstyle T}} V^{-1} y\right)^{\mathrm{\scriptscriptstyle T}} X^{\mathrm{\scriptscriptstyle T}} Q^{-1}\left(y-m_0\right) \\ & =y^{\mathrm{\scriptscriptstyle T}} V^{-1}\left(y-X m_0\right)-\left(X m_0+X M_0 X^{\mathrm{\scriptscriptstyle T}} V^{-1} y\right)^{\mathrm{\scriptscriptstyle T}} Q^{-1}\left(y-X m_0\right)\\ & =y^{\mathrm{\scriptscriptstyle T}} V^{-1}\left(y-X m_0\right)-\left(Q^{-1} X m_0+Q^{-1}\left(X M_0 X^{\mathrm{\scriptscriptstyle T}}\right)V^{-1} y\right)\left(y-X m_0\right) \\ & =y^{\mathrm{\scriptscriptstyle T}} V^{-1}\left(y-X m_0\right)-\left(Q^{-1} X m_0+Q^{-1}(Q-V) V^{-1} y\right)^{\mathrm{\scriptscriptstyle T}}(y-X m_0) \\ & =y^{\mathrm{\scriptscriptstyle T}} V^{-1}\left(y-X m_0\right)-\left(Q^{-1} X m_0+V^{-1} y- Q^{-1} y \right)^{\mathrm{\scriptscriptstyle T}}\left(y-X m_0\right) \\ & =y^{\mathrm{\scriptscriptstyle T}} V^{-1}\left(y-X m_0\right)-\left(V^{-1} y+Q^{-1}\left(X m_0-y\right)\right)^{\mathrm{\scriptscriptstyle T}}\left(y-X m_0\right) \\ & =y^{\mathrm{\scriptscriptstyle T}} V^{-1}\left(y-X m_0\right)-y^{\mathrm{\scriptscriptstyle T}} V^{-1}\left(y-X m_0\right) +\left(y-X m_0\right)^{\mathrm{\scriptscriptstyle T}} Q^{-1}\left(y-X m_0\right) \\ & =\left(y-X m_0\right)^{\mathrm{\scriptscriptstyle T}} Q^{-1}\left(y-X m_0\right) \\ \end{align}$$

So, we get the following updating form of the posterior distribution from Bayesian linear regression

$$\begin{aligned} P\left(\beta, \sigma^{2} \mid y\right) = NIG\left(\beta, \sigma^{2} \mid \tilde{m}_1, \tilde{M}_1, a_{1}, b_{1}\right) \; \end{aligned}$$ where

$$\begin{aligned} \tilde{m}_1 & =M_1 m_1=m_0+M_0 X^{\mathrm{\scriptscriptstyle T}} Q^{-1}\left(y-X m_0\right) \\ \tilde{M}_1 & =M_1=M_0-M_0 X^{\mathrm{\scriptscriptstyle T}} Q^{-1} X M_0 \\ a_1 & =a_0+\frac{p}{2} \\ b_1 & =b_0+\frac{1}{2}\left(y-X m_0\right)^{\mathrm{\scriptscriptstyle T}} Q^{-1}\left(y-X m_0\right) \\ Q & =V+X M_0 X^{\mathrm{\scriptscriptstyle T}} \end{aligned}$$

1.4.2 Method 2: Distribution theory

Previously, we got the Bayesian Linear Regression Updater using Sherman-Woodbury-Morrison identity. Here, we will derive the results without resorting to it. Recall that the model is given by

$$\begin{align} & y=X \beta+\epsilon \ , \ \epsilon \sim N\left(0, \sigma^2 V\right) \\ & \beta=m_0+\omega \ , \ \omega \sim N\left(0, \sigma^2 M_0\right) \\ & \sigma^2 \sim I G\left(a_0, b_0\right) \end{align}$$

This corresponds to the posterior distribution

$$\begin{align} P\left(\beta, \sigma^2 \mid y\right) \propto I G\left(\sigma^2 \mid a_0, b_0\right) & \times N\left(\beta \mid m_0, \sigma^2 M_0\right) \times N\left(y \mid X \beta, \sigma^2 V\right) \end{align}$$

We will derive $P\left(\sigma^2 \mid y\right)$ and $P\left(\beta \mid \sigma^2, y\right)$ in a form that will reflect updates from the prior to the posterior. Integrating out $\beta$ from the model is equivalent to substituting $\beta$ from its prior model. Thus, $P\left(y \mid \sigma^2\right)$ is derived simply from $y =X \beta+\epsilon =X\left(m_0+\omega\right)+\epsilon =X m_0 + X \omega + \epsilon =X m_0+ \eta$

where $\eta = X \omega + \epsilon \sim N\left(0, \sigma^2Q\right)$ and $Q=X M_0 X^{\mathrm{\scriptscriptstyle T}}+V \, .$

Therefore,

$$\begin{align} y \mid \sigma^2 \sim N\left(X m_0, \sigma^2 Q\right)\\ \end{align}$$

The posterior distribution is given by $$\begin{align} P\left(\sigma^2 \mid y\right) & \propto IG \left(\sigma^2 \mid a_1, b_1\right) \end{align}$$

where $$\begin{align} & a_1 = a_0 + \frac{p}{2} \\ & b_1 = b_0 + \frac{1}{2} (y-Xm_0)^{\mathrm{\scriptscriptstyle T}} Q^{-1} \left(y-Xm_0\right) \end{align}$$

$$\begin{align} P\left(\sigma^2 \mid y\right) & \propto P\left(\sigma^2\right) P\left(y \mid \sigma^2\right) \\ & \propto\left(\sigma^2 \mid a_0, b_0\right) \times N\left(y \mid X m_0, \sigma^2 Q\right) \\ & \propto\left(\frac{1}{\sigma^2}\right)^{a_0+1} e^{-\frac{b_0} {\sigma^2} \times\left(\frac{1}{\sigma^2}\right)^{\frac{n}{2}} e^{-\frac{1}{2 \sigma^2}}\left(y-Xm_0\right)^{\mathrm{\scriptscriptstyle T}} Q^{-1}\left(y-Xm_0\right)} \\ & \propto\left(\frac{1}{\sigma^2}\right)^{a_0+\frac{p}{2}+1} e^{-\frac{1}{\sigma^2}\left(b_0+\frac{1}{2}\left(y-Xm_0\right)^{\mathrm{\scriptscriptstyle T}} Q^{-1}\left(y-Xm_0\right)\right)} \\ & \propto IG \left(\sigma^2 \mid a_1, b_1\right) \end{align}$$

where $$\begin{align} & a_1 = a_0 + \frac{p}{2} \\ & b_1 = b_0 + \frac{1}{2} (y-Xm_0)^{\mathrm{\scriptscriptstyle T}} Q^{-1} \left(y-Xm_0\right) \end{align}$$

Next, we turn to $P\left(\beta \mid \sigma^2, y\right)$. Note that $$\left(\begin{array}{l} y \\ \beta \end{array}\right) \mid \sigma^2 \sim N\left(\left(\begin{array}{l} Xm_0 \\ m_0 \end{array}\right), \quad \sigma^2 \left(\begin{array}{cc} Q & X M_0 \\ M_0 X^{\mathrm{\scriptscriptstyle T}} & M_0 \end{array}\right)\right) \;$$

We have used the facts

$$\begin{aligned} & \operatorname{E}\left[y \mid \sigma^2\right] = Xm_0 \, , \ \operatorname{Var}\left(y \mid \sigma^2\right)=\sigma^2 Q \, ; \\ & \operatorname{E}\left[\beta \mid \sigma^2\right] = m_0 \, , \ \operatorname{Var}\left(\beta \mid \sigma^2\right)=\sigma^2 M_0 \, ; \end{aligned}$$

$$\begin{aligned} \operatorname{Cov}\left(y, \beta \mid \sigma^2\right) &= \operatorname{Cov}\left(X \beta+\epsilon, \beta \mid \sigma^2\right) \\ & =\operatorname{Cov}\left(X\left(m_0+\omega\right)+\epsilon, m_0+\omega \mid \sigma^2\right) \\ & \quad \left( \text {Since } m_0 \text { is constant and } \operatorname{Cov}(\omega, \epsilon)=0 \right) \\ & =\operatorname{Cov}\left(X \omega, \omega \mid \sigma^2\right) \\ & =\sigma^2 X M_0 \end{aligned}$$

From the expression of a conditional distribution derived from a multivariate Gaussian, we obtain $$\beta \mid \sigma^2, y \sim N\left(\tilde{m}_1, \sigma^2 \tilde{M}_1\right)$$

where $$\begin{align} & \tilde{m}_1=\operatorname{E}\left[\beta \mid \sigma^2, y\right]=m_0+M_0 X^{\mathrm{\scriptscriptstyle T}} Q^{-1}\left(y-X{m_0}\right) \\ & \tilde{M}_1=M_0-M_0 X^{\mathrm{\scriptscriptstyle T}} Q^{-1} X M_0 \\ \end{align}$$

Note:

$$\begin{align} & \left(\begin{array}{l} X_1 \\ X_2 \end{array}\right) \sim N\left(\left(\begin{array}{l} \mu_1 \\ \mu_2 \end{array}\right),\left(\begin{array}{ll} \Sigma_{11} & \Sigma_{12} \\ \Sigma_{21} & \Sigma_{22} \end{array}\right)\right) \text { with } \Sigma_{21} = \Sigma_{12}^{\mathrm{\scriptscriptstyle T}} \end{align}$$

$$\begin{align} & \Rightarrow X_2 \mid X_1 \sim N\left(\mu_{2 \cdot 1}, \Sigma_{2 \cdot 1}\right) \end{align}$$

where $\mu_{2 \cdot 1}= \mu_2+\Sigma_{21} \Sigma_{11}^{-1}\left(X_1-\mu_1\right) \text { and } \Sigma_{2 \cdot 1}=\Sigma_{22}-\Sigma_{21} \Sigma_{11}^{-1} \Sigma_{12} \, .$

1.5 Bayesian prediction

Assume $V=I_{n}$. Let $\tilde{y}$ denote an $\tilde{n}\times 1$ vector of outcomes. $\tilde{X}$ is corresponding predictors. We seek to predict $\tilde{y}$ based upon $y$

$$\begin{align} P(\tilde{y} \mid y) &= t_{2a_1}\left(\tilde{X} M_1 m_1, \frac{b_1}{a_1}\left(I_{\tilde{n}} + \tilde{X} M_{1} \tilde{X}^{\mathrm{\scriptscriptstyle T}}\right)\right) \; \end{align}$$

$$\begin{align} P\left(\beta, \sigma^{2}, \tilde{y} \mid y\right) &=P\left(\beta, \sigma^{2} \mid y\right) \cdot P\left(\tilde{y} \mid \beta, \sigma^{2}, y\right) \\ &= P\left(\beta, \sigma^{2}\right) \cdot P\left(y \mid \beta, \sigma^{2}\right) \cdot P\left(\tilde{y} \mid \beta, \sigma^{2}, y\right) \\ &= NIG \left(\beta, \sigma^{2} \mid m_{0}, M_{0}, a_{0}, b_{0}\right) \cdot N\left(y \mid X \beta, \sigma^{2} I_{n}\right) \cdot N\left(\tilde{y} \mid \tilde{X} \beta, \sigma^{2} I_{\tilde{n}}\right) \\ &= NIG \left(\beta, \sigma^{2} \mid M_{1} m_{1}, M_{1}, a_{1}, b_{1}\right) \cdot N\left(\tilde{y} \mid \tilde{X} \beta, \sigma^{2} I_{\tilde{n}}\right) \\ &= IG\left(\sigma^{2} \mid a_{1}, b_{1}\right) \cdot N\left(\beta \mid M_{1} m_{1}, \sigma^{2} M_{1} \right) \cdot N\left(\tilde{y} \mid \tilde{X} \beta, \sigma^{2} I_{\tilde{n}} \right) \; \end{align}$$

Then we can calculate posterior predictive density $P\left(\tilde{y} \mid y\right)$ from $P\left(\beta, \sigma^{2}, \tilde{y} \mid y\right)$

$$\begin{align} P\left(\tilde{y} \mid y\right) &=\iint P\left(\beta, \sigma^{2}, \tilde{y} \mid y\right) \ d\beta \ d\sigma^{2} \\ &=\iint IG\left(\sigma^{2} \mid a_{1}, b_{1}\right) \cdot N\left(\beta \mid M_{1} m_{1}, \sigma^{2} M_{1} \right) \cdot N\left(\tilde{y} \mid \tilde{X} \beta, \sigma^{2} I_{\tilde{n}}\right) \ d\beta \ d\sigma^{2} \\ &=\int IG\left(\sigma^{2} \mid a_{1}, b_{1}\right) \int N\left(\beta \mid M_{1} m_{1}, \sigma^{2} M_{1} \right) \cdot N\left(\tilde{y} \mid \tilde{X} \beta, \sigma^{2} I_{\tilde{n}}\right) \ d\beta \ d\sigma^{2} \\ \end{align}$$

As for $\int N\left(\beta \mid M_{1} m_{1}, \sigma^{2} M_{1}\right) \cdot N\left(\tilde{y} \mid \tilde{X} \beta, \sigma^{2} I_{\tilde{n}}\right) \ d\beta$, we provide an easy way to derive it avoiding any integration at all. Note that we can write the above model as

$$\begin{align} \tilde{y} &= \tilde{X} \beta + \tilde{\epsilon}, \text{ where } \tilde{\epsilon} \sim N\left(0,\sigma^2 I_{\tilde{n}}\right) \\ \beta &= M_{1} m_{1} + \epsilon_{\beta \mid y}, \text{ where } \epsilon_{\beta \mid y} \sim N\left(0,\sigma^2M_{1}\right) \end{align}$$

where $\tilde{\epsilon}$ and $\epsilon_{\beta \mid y}$ are independent of each other. It then follows that

$$\begin{align} \tilde{y} &= \tilde{X} M_{1} m_{1} + \tilde{X} \epsilon_{\beta \mid y} + \tilde{\epsilon} \sim N\left(\tilde{X} M_{1} m_{1}, \sigma^2\left(I_{\tilde{n}} + \tilde{X} M_{1} \tilde{X}^{\mathrm{\scriptscriptstyle T}}\right)\right) \end{align}$$

As a result

$$\begin{align} P\left(\tilde{y} \mid y\right) &=\int IG\left(\sigma^{2} \mid a_{1}, b_{1}\right) \cdot N\left(\tilde{X} M_{1} m_{1}, \sigma^2\left(I_{\tilde{n}} + \tilde{X} M_{1} \tilde{X}^{\mathrm{\scriptscriptstyle T}}\right)\right) \ d\sigma^{2} \\ &= t_{2a_1}\left(\tilde{X} M_1 m_1, \frac{b_1}{a_1}\left(I_{\tilde{n}} + \tilde{X} M_{1} \tilde{X}^{\mathrm{\scriptscriptstyle T}}\right)\right) \; \end{align}$$

1.6 Sampling from the posterior distribution

We can get the joint posterior density $P\left(\beta, \sigma^{2}, \tilde{y} \mid y\right)$ by sampling process

1. Draw $\hat{\sigma}_{(i)}^{2}$ from $I G\left(a_{1}, b_{1}\right)$

2. Draw $\hat{\beta}_{(i)}$ from $N\left(M_{1} m_{1}, \hat{\sigma}_{(i)}^{2} M_{1}\right)$

3. Draw $\tilde{y}_{(i)}$ from $N\left(\tilde{X} \hat{\beta}_{(i)}, \hat{\sigma}_{(i)}^{2}I_{\tilde{n}}\right)$