如果移动端无法正确渲染某些公式,建议转PC站点,也是做了自适应处理的。点击这里跳转

BMSI系列的第七篇,介绍叶斯框架下的线性回归

手敲$\LaTeX$难免出现纰漏,有任何疑似错误或者不清楚的地方请直接在下方评论区留言(评论区只开在【Series】下),谢谢各位读者。

线性回归是统计中极其基础而十分重要的部分,相信读者对其不会陌生。本节则探讨贝叶斯框架下的线性回归的形式,并和传统的线性回归进行对比加深理解。

传统的线性回归的结论

我们使用记号$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)$$