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)\]