BMSI系列的第七篇，介绍叶斯框架下的线性回归
手敲\(\LaTeX\)难免出现纰漏，有任何疑似错误或者不清楚的地方请直接在下方评论区留言，谢谢各位读者。
线性回归是统计中极其基础而十分重要的部分，相信读者对其不会陌生。本节则探讨贝叶斯框架下的线性回归的形式，并和传统的线性回归进行对比加深理解。
传统的线性回归的结论
我们使用记号\(y\lvert X\sim\mathcal{N}(X\beta,\sigma^2 I)\)。根据数理统计的知识不难得参数\(\beta\)的极大似然估计为
\[\widehat{\beta}=\widehat{\beta}_{\mathrm{MLE}}=(X^\intercal X)^{-1}X^\intercal y\]
且其分布为
\[\widehat{\beta}\sim t_{n-k}(\beta, s^2(X^\intercal X)^{-1})\]
其中\(s^2=\mathrm{SSE}/(n-k),\,\mathrm{SSE}=(Y-X\widehat{\beta})^\intercal(Y-X\widehat{\beta})\)
而\(s^2\)的分布为
\[\frac{(n-k)s^2}{\sigma^2}\sim\chi_{n-k}^2\ \Rightarrow\ \frac{1}{s^2}\sim\mathrm{Inv-}\chi^2(n-k,\sigma_0^{-2})\]
模型假设
我们假设，自变量\(X\sim p(X\lvert\Phi)\)，因变量\(Y\lvert X\sim p(y\lvert X,\theta)\)，进一步地，我们认为\(p(\Phi,\theta)=p(\Phi)p(\theta)\)即二者独立。那么模型参数的后验分布
\[p(\Phi,\theta\lvert X,y)=\frac{p(X\lvert\Phi)p(\Phi)p(y\lvert X,\theta)p(\theta)}{p(X)p(y\lvert X)}=p(\Phi\lvert X)p(\theta\lvert X,y)\]
可以看到，如果我们只关心\(\theta\)的话，我们只需要关注\(p(\theta\lvert X,y)\)，这和我们传统的线性回归是一致的。而
\[p(\theta\lvert X,y)\propto p(y\lvert X,\theta)p(\theta)\]
为方便记，下文中将\(X,y\)合并为\(y\)表示数据信息。
无信息先验
推导
此时，\(p(\theta)=p(\beta,\sigma^2)\propto\sigma^{-2}\)。现在我们求取其后验。类似于第二节中多参数模型的分析，我们先\(p(\beta,\sigma^2\lvert y)=p(\beta\lvert\sigma^2,y)p(\sigma^2\lvert y)\)，求出后两个后再求边缘密度\(p(\beta\lvert y)\)。先看第一项，有
\[p(\beta|\sigma^2,y)\propto p(y|\beta,\sigma^2)p(\beta)\propto \exp\{-\frac{(Y-X\beta)^\intercal(Y-X\beta)}{2\sigma^2} \}\]
可以看到指数部分是\(\beta\)的二次型，因此设\(\beta\sim\mathcal{N}(\widehat{\beta},V_\beta\sigma^2)\)，于是有
\[\frac{(Y-X\beta)^\intercal (Y-X\beta)}{\sigma^2} = \frac{(\beta-\widehat{\beta})^\intercal V_{\beta}^{-1}(\beta-\widehat{\beta})}{\sigma^2}\]
对比相关项我们有
\[X^\intercal Y=V_{\beta}^{-1}\widehat{\beta},\quad X^\intercal X=V_{\beta}^{-1}\]
也就是
\[\widehat{\beta}=(X^\intercal X)^{-1}X^\intercal Y,\quad V_{\beta}=(X^\intercal X)^{-1}\]
接下来求\(p(\sigma^2\lvert y)\)，我们有
\[ p(\sigma^2|y)=\frac{p(\beta, \sigma^2|y)}{p(\beta|\sigma^2,y)}=\frac{\sigma^{-2}\sigma^{-n}\exp\{-\frac{(Y-X\beta)^\intercal (Y-X\beta)}{2\sigma^2} \}}{\sqrt{|(X^\intercal X)^{-1}\sigma^2|}} = (\sigma^2)^{-(n-k)/2-1}\exp\big\{-\frac{(n-k)s^2}{2\sigma^2}\big\} \]
也就是\(\sigma^2\sim\mathrm{Inv-}\chi^2(n-k,s^2)\)
由此，进一步得到\(\beta\lvert y\sim t_{n-k}(\widehat{\beta},s^2V_\beta)\)（参见多参数模型）。但事实上我们很少使用这一结论，因为在抽样的过程中，我们完全可以先抽\(\sigma\)，再在此基础上抽\(\beta\)
对比
可以看到二者的形式极为相似，让我们更加细致地查看二者的差别
在传统的回归中，我们有真实、确定但未知的参数\(\beta_0,\sigma_0^2\)，参数的不确定性是由于使用的样本不能完全代表整体，也就是抽样时所引入的。此时
\[1/s^2\sim\mathrm{Inv-}\chi^2(n-k,1/\sigma_0^2)\] \[\widehat{\beta}\sim t_{n-k}(\beta_0, s^2(X^\intercal X)^{-1})\]
而在无信息先验的贝叶斯线性回归中，参数的不确定性由后验分布来解释
\[\sigma^2\lvert y\sim\mathrm{Inv-}\chi^2(n-k,s^2)\] \[\beta\lvert y\sim t_{n-k}(\widehat{\beta}, s^2(X^\intercal X)^{-1})\]
共轭先验
推导
我们取
\[\beta\lvert\sigma^2\sim\mathcal{N}(m_0,\sigma^2C_0),\ \sigma^2\sim\mathrm{Inv-}\chi^2(v_0,s_0^2)\]
那么后验分布为
\[\beta\lvert\sigma^2,y\sim\mathcal{N}(m_n,\sigma^2C_n),\ \sigma^2\lvert y\sim\mathrm{Inv-}\chi^2(v_n,s_n^2)\]
其中
\[\begin{align} m_n&=m_0+C_0X^\intercal(XC_0X^\intercal+I)^{-1}(y-Xm_0)\\ C_n&=C_0-C_0X^\intercal(XC_0X^\intercal+I)^{-1}XC_0\\ v_n&=v_0+n\\ v_ns_n^2&=v_0s_0^2+(y-Xm_0)^\intercal(XC_0X^\intercal+I)^{-1}(y-Xm_0)\\ \end{align}\]
在R中，arm包里的bayesglm提供了共轭先验下的贝叶斯线性回归的一个实现，可以通过设定priro.mean等信息调整先验，具体的请使用??arm::bayesglm查看文档。
数值计算的技巧
可以看到许多地方都出现了\((X^\intercal X)^{-1}\)。我们知道矩阵求逆向来都是计算的瓶颈。一个优化的技巧是使用QR分解。此时
\[V_\beta=(X^\intercal X)^{-1}=(R^\intercal R)^{-1}=R^{-1}(R^\intercal)^{-1}\] \[R\widehat{\beta}=R(X^\intercal X)^{-1}X\beta=RR^{-1}(R^\intercal)^{-1}R^\intercal Q^\intercal y=Q^\intercal y\]
注意到\(R\)是上三角阵，因此我们无需求\(R\)的逆。
特殊\(\beta\)的先验
我们可以选取形式为\(\beta\lvert\sigma^2\sim\mathcal{N}(b,\sigma^2 B)\)的先验。其中有一些可以考虑的特殊情况
- \(b=0\)
- \(B\) 为对角阵
- \(B=gI\)
- \(B=g(X^\intercal X)^{-1}\)
其中，第四种情况也称为Zellner's g-prior。它实现了先验和数据（似然）的一个折衷。对于\(g\)的选取有以下五种考虑
- \(g=1\)意味着模型给了先验和数据（总体）相同的权重
- \(g=n\)意味着先验的权重相当于一个数据点的权重
- \(g=\infty\)意味着无信息先验
- \(g=\widehat{g}_{\mathrm{EG}}=\arg\max_gp(y\lvert g)\)，即Empirical Bayes estimate。其中
\[p(y\lvert g)=\frac{\Gamma((n-1)/2)}{\sqrt{n\pi^{n+1}}}\lvert\lvert y-\bar{y}\rvert\rvert^{-(n-1)}\frac{(1+g)^{(n-1-k)/2}}{(1+g(1-R^2))^{(n-1)/2}}\] 其中\(R^2\)为回归的决定系数
- 也给\(g\)一个先验，然后做一个完整的贝叶斯分析
在R中，BMS包里的zlm提供了Zellner's g-prior下的贝叶斯线性回归的一个实现，可以通过设定g来调节结果，具体的请使用??BMS::zlm查看文档。
不同的方差结构
我们之前的讨论中都是假设因变量的方差是独立且相等的，也就是说协方差是单位阵（乘一个标量）。而现在我们放宽这一要求，即\(y\lvert \beta,\Sigma_y\sim\mathcal{N}(X\beta,\Sigma_y)\)。当然，我们是要求\(\Sigma_y\)是一个对称正定阵。
确定的方差结构
实际使用中，我们通常将\(\Sigma_y\)设定为\(\sigma^2Q_y\)，而\(Q_y\)是完全给定的，\(\sigma^2\)则并不知道。也就是我们给定了方差的结构。一个平凡的例子自然是\(Q_y=I\)，也就是之前讨论过的。
在这样的框架下，我们有
\[Q_y^{-1/2}y\lvert\beta,\sigma^2\sim\mathcal{N}(Q_y^{-1/2}X\beta,\sigma^2I)\]
也就化归为之前讨论过的情形。使用无信息先验\(p(\beta,\sigma^2)\propto\sigma^{-2}\)，有
\[\begin{align} \widehat{\beta}&=(X^\intercal Q_y^{-1}X)^{-1}X^\intercal Q_y^{-1}y\\ V_\beta&=(X^\intercal Q_y^{-1}X)^{-1}\\ s^2=(y-X\widehat{\beta})^\intercal Q_y^{-1}(y-X\widehat{\beta})/(n-k) \end{align}\]
计算时，由于\(Q_y\)是对称正定阵，我们可以使用Cholesky分解或奇异值分解求得\(Q_y^{1/2}\)以及\(Q_y^{-1/2}\)
几个常见的\(Q_y\)有
- \(\mathrm{diag}(1/w_1,1/w_2,\ldots,1/w_n)\)。这通常会发生在\(y_i\)是\(w_i\)个观测的平均时
- 等相关阵，如\(Q_y(\rho)=\left(\begin{matrix}1&\rho&\rho\\\rho&1&\rho\\\rho&\rho&1\end{matrix}\right)\)
- AR(1)，如\(Q_y(\rho)=\left(\begin{matrix}1&\rho&\rho^2\\\rho&1&\rho\\\rho^2&\rho&1\end{matrix}\right)\)
注意到，后两个矩阵其实依赖于新的参数\(\rho\)，我们将其记为\(\Phi\)。同样的，我们使用无信息先验\(p(\beta,\sigma^2)\propto\sigma^{-2}\)，有
\[p(\Phi,y)=\frac{p(\beta,\sigma^2,\Phi\lvert y)}{p(\beta,\sigma^2\lvert\Phi,y)}\propto\frac{p(\Phi)\mathcal{N}(y;X\beta,\sigma^2Q_y)}{\mathrm{Inv-}\chi^2(\sigma^2;n-k,s^2)\mathcal{N}(\beta;\widehat{\beta},V_\beta\sigma^2)}\propto p(\Phi)\lvert V_\beta\rvert^{1/2}(s^2)^{-(n-k)/2}\]
不定的方差结构
我们也可以不指定方差的结构，只要满足对称正定性就可以。这里我们讨论两个先验
\(p(\beta\lvert\Sigma_y)\propto 1\)
此时
\[\beta\lvert\Sigma_y,y\sim\mathcal{N}((X^\intercal\Sigma_y^{-1}X)^{-1}X^\intercal y,(X^\intercal\Sigma_y^{-1}X)^{-1})\]
\[p(\Sigma_y\lvert y)\propto p(\Sigma_y)\lvert(X^\intercal\Sigma_y^{-1}X)\rvert^{-1/2}\exp\left(-\frac{1}{2}(y-X\widehat{\beta})^\intercal\Sigma_y^{-1}(y-X\widehat{\beta})\right)\]
此时，我们熟悉的\(\Sigma_y\sim\mathrm{Inv-Wishart}_v(S^{-1})\)是其共轭先验。
\(\beta\lvert\Sigma_y\sim\mathcal{N}(\beta_0,\Sigma_\beta)\)
此时，\(\beta\lvert\Sigma_y,y\sim\mathcal{N}(\mu,\Lambda)\)
其中
\[\Lambda=(\Sigma_\beta^{-1}+X^\intercal\Sigma_y^{-1}X)^{-1}\] \[\mu=\Lambda(\Sigma_\beta^{-1}\beta_0+X^\intercal\Sigma_y^{-1}X)\]