线性回归系数的几个性质

摘要
线性回归问题的描述
单变量线性回归系数的公式
无偏估计
线性回归系数的方差
其余的几个性质
- 残差项之和为0
- 线性拟合直线总会经过 (xˉ,yˉ)(\bar{x}, \bar{y})(xˉ,yˉ) 这个点
- 在 xix_ixi 权重下，残差和为0，即∑i=1nxiei=0\displaystyle \sum_{i = 1}^n x_i e_i = 0i=1∑nxiei=0
- 在 yi^\hat{y_i}yi^ 权重下，残差和为0，即∑i=1nyi^ei=0\displaystyle \sum_{i = 1}^n \hat{y_i} e_i = 0i=1∑nyi^ei=0
模拟
多元线性回归的系数也是无偏估计
参考文献

摘要

我们知道一元和多元线性回归系数都有解析解，本文将简要介绍总结线性回归系数的几个常见的性质。

线性回归问题的描述

我们回忆一下，单变量线性回归问题是指，给定了 nnn 个观察量 (x1,y1),(x2,y2),⋯,(xn,yn)(x_1, \, y_1), \, (x_2, \, y_2), \cdots, (x_n, y_n)(x1,y1),(x2,y2),⋯,(xn,yn)。我们希望用一个线性的关系 y=β1x+β0+ϵy = \beta_1 x + \beta_0 + \epsilony=β1x+β0+ϵ 来描述这些观察量的规律。

这里，我们把方程
y=β1x+β0+ϵy = \beta_1 x + \beta_0 + \epsilony=β1x+β0+ϵ

称为 总体回归模型 (population regression model)。

而当给定的 nnn 个观察量 (x1,y1),(x2,y2),⋯,(xn,yn)(x_1, \, y_1), \, (x_2, \, y_2), \cdots, (x_n, y_n)(x1,y1),(x2,y2),⋯,(xn,yn)，我们称方程
yi=β1xi+β0+ϵiy_i = \beta_1 x_i + \beta_0 + \epsilon_iyi=β1xi+β0+ϵi
为 样本回归模型（sample regression model）。
其中 ϵi\epsilon_iϵi 为误差项 (error)，每个 ϵi\epsilon_iϵi 均是一个随机变量，独立且均服从一个均值为0, 方差为 σ2\sigma^2σ2 的概率分布。其中 σ2\sigma^2σ2 为误差项的方差，在我们所考虑的线性模型中，我们已知 σ2\sigma^2σ2 是固定的，但是我们不知道其具体的数值。

所以我们去做线性回归“拟合”模型的参数 β1\beta_1β1 和 β0\beta_0β0 时，实际上是根据 (x1,y1),(x2,y2),⋯,(xn,yn)(x_1, \, y_1), \, (x_2, \, y_2), \cdots, (x_n, y_n)(x1,y1),(x2,y2),⋯,(xn,yn) 去做点估计 (point estimator)。

这里值得注意的是，我们把 x1,x2,⋯,xnx_1, \, x_2, \, \cdots, \, x_nx1,x2,⋯,xn 当作给定的值，即可以认为是常量，而每一个 yi,i=1,2,⋯ny_i, \, i = 1, 2, \, \cdots \, nyi,i=1,2,⋯n 均是一个随机变量。

所以， yiy_iyi 作为随机变量，可以用下图来体现。

上图中，xxx 轴上的绿色的圆点表示固定的 xix_ixi ，而橙色的直线表示真实的线性关系，而蓝色的钟型曲线代表随机变量 yiy_iyi 服从的概率密度函数。

单变量线性回归系数的公式

我们回顾单变量线性回归问题的公式，有
{β1^=∑(xi−xˉ)(yi−yˉ)∑(xi−xˉ)2β0^=yˉ−β1^xˉ\begin{cases} & \hat{\beta_1} = \frac{\sum (x_i - \bar{x} ) (y_i - \bar{y})}{\sum (x_i - \bar{x} )^2 } \\ & \hat{\beta_0} = \bar{y} - \hat{\beta_1} \bar{x} \\ \end{cases} {β1^=∑(xi−xˉ)2∑(xi−xˉ)(yi−yˉ)β0^=yˉ−β1^xˉ

即我们对参数 β1,β0\beta_1, \, \beta_0β1,β0 的估计是 β1^=∑(xi−xˉ)(yi−yˉ)∑(xi−xˉ)2\displaystyle \hat{\beta_1} = \frac{\sum (x_i - \bar{x} ) (y_i - \bar{y})}{\sum (x_i - \bar{x} )^2 }β1^=∑(xi−xˉ)2∑(xi−xˉ)(yi−yˉ)，对参数 β0\beta_0β0 的估计是 β0^=yˉ−β1^xˉ\displaystyle \hat{\beta_0} = \bar{y} - \hat{\beta_1} \bar{x}β0^=yˉ−β1^xˉ。

为方便起见，我们记

Sxx=∑i=1n(xi−xˉ)2\displaystyle S_{xx} = \sum_{i = 1}^n (x_i - \bar{x})^2Sxx=i=1∑n(xi−xˉ)2,
Sxy=∑i=1n(xi−xˉ)(yi−yˉ)\displaystyle S_{xy} = \sum_{i = 1}^n (x_i - \bar{x}) (y_i - \bar{y})Sxy=i=1∑n(xi−xˉ)(yi−yˉ)。

我们可以进一步简化 SxyS_{xy}Sxy 为 Sxy=∑i=1n(xi−xˉ)yi\displaystyle S_{xy} = \sum_{i = 1}^n (x_i - \bar{x}) y_iSxy=i=1∑n(xi−xˉ)yi。

那么我们可以将 β1^\hat{\beta_1}β1^ 与 β0^\hat{\beta_0}β0^ 写成：
β1^=∑i=1n(xi−xˉ)Sxx⋅yi\displaystyle \hat{\beta_1} =\sum_{i = 1}^n \frac{ (x_i - \bar{x} ) }{S_{xx} } \cdot y_iβ1^=i=1∑nSxx(xi−xˉ)⋅yi,
β0^=yˉ−β1^xˉ\displaystyle \hat{\beta_0} = \bar{y} - \hat{\beta_1} \bar{x}β0^=yˉ−β1^xˉ。

回忆起之前我们提到的，每一个 yi,i=1,2,⋯ny_i, \, i = 1, 2, \, \cdots \, nyi,i=1,2,⋯n 均是一个随机变量，这里我们分别对 β1\beta_1β1 和 β0\beta_0β0 的估计 β1^\hat{\beta_1}β1^ 与 β0^\hat{\beta_0}β0^ 就是 yi,i=1,2,⋯ny_i, i = 1, \, 2, \, \cdots \, nyi,i=1,2,⋯n 的函数。

无偏估计

首先，我们证明上述估计 β1^\hat{\beta_1}β1^ 与 β0^\hat{\beta_0}β0^ 是无偏估计 (unbiased)。

要证明我们的估计是无偏估计，我们须要证明我们估计的期望恒等于所估计的参数，即我们须要证明：
{E[β1^]=β1E[β0^]=β0\begin{cases} & \mathbb{E}[ \hat{\beta_1} ] = \beta_1 \\ \\ & \mathbb{E}[ \hat{\beta_0} ] = \beta_0 \\ \end{cases} ⎩⎨⎧E[β1^]=β1E[β0^]=β0

证明过程十分直接，我们直接将上一节的表达式代入。

E[β1^]=E[∑i=1n(xi−xˉ)Sxx⋅yi]\displaystyle \mathbb{E}[ \hat{\beta_1} ] =\mathbb{E} \Big[ \sum_{i = 1}^n \frac{ (x_i - \bar{x} ) }{S_{xx} } \cdot y_i \Big]E[β1^]=E[i=1∑nSxx(xi−xˉ)⋅yi]，注意到 yi=β1xi+β0+ϵiy_i = \beta_1 x_i + \beta_0 + \epsilon_iyi=β1xi+β0+ϵi，我们有 E(yi)=β1xi+β0\displaystyle \mathbb{E}(y_i) = \beta_1 x_i + \beta_0E(yi)=β1xi+β0。这里我们用到了 E[ϵi]=0\displaystyle \mathbb{E} [ \epsilon_i ] = 0E[ϵi]=0。

代入，我们有

E[β1^]=E[∑i=1n(xi−xˉ)Sxx⋅yi]=∑i=1n(xi−xˉ)SxxE[yi]=∑i=1n(xi−xˉ)Sxx⋅(β1xi+β0)=β1⋅[∑i=1n(xi−xˉ)Sxx⋅xi]+β0⋅∑i=1n(xi−xˉ)Sxx=β1\begin{aligned} \displaystyle \mathbb{E} [\hat{\beta_1}] &= \mathbb{E} \Big[ \sum_{i = 1}^n \frac{ (x_i - \bar{x} ) }{S_{xx} } \cdot y_i \Big] \\ &= \sum_{i = 1}^n \frac{ (x_i - \bar{x} ) }{S_{xx}} \mathbb{E} [y_i] \\ &= \sum_{i = 1}^n \frac{ (x_i - \bar{x} ) }{S_{xx} } \cdot (\beta_1 x_i + \beta_0) \\ &= \beta_1 \cdot \Big[ \sum_{i = 1}^n \frac{ (x_i - \bar{x} ) }{S_{xx}} \cdot x_i \Big] + \beta_0 \cdot \sum_{i = 1}^n \frac{ (x_i - \bar{x} ) }{S_{xx}} \\ &= \beta_1 \end{aligned}E[β1^]=E[i=1∑nSxx(xi−xˉ)⋅yi]=i=1∑nSxx(xi−xˉ)E[yi]=i=1∑nSxx(xi−xˉ)⋅(β1xi+β0)=β1⋅[i=1∑nSxx(xi−xˉ)⋅xi]+β0⋅i=1∑nSxx(xi−xˉ)=β1。

注意，上式中，我们用到了 ∑i=1n(xi−xˉ)Sxx=0\displaystyle \sum_{i = 1}^n \frac{ (x_i - \bar{x} ) }{S_{xx}} = 0i=1∑nSxx(xi−xˉ)=0，以及 ∑i=1n(xi−xˉ)⋅xi=Sxx\displaystyle \sum_{i = 1}^n (x_i - \bar{x}) \cdot x_i = S_{xx}i=1∑n(xi−xˉ)⋅xi=Sxx。

同样的，我们有，
E(β0^)=E[yˉ−xˉ⋅β1^]=1n∑i=1n(β0+β1⋅xi)−xˉ⋅β1=β0\begin{aligned} \displaystyle \mathbb{E} (\hat{\beta_0}) &= \mathbb{E} \big[ \bar{y} - \bar{x} \cdot \hat{\beta_1} \big] \\ &= \frac{1}{n} \sum_{i = 1}^n (\beta_0 + \beta_1 \cdot x_i) - \bar{x} \cdot \beta_1 \\ &= \beta_0 \end{aligned}E(β0^)=E[yˉ−xˉ⋅β1^]=n1i=1∑n(β0+β1⋅xi)−xˉ⋅β1=β0。

所以，β0^\hat{\beta_0}β0^ 也是 β0\beta_0β0 的无偏估计。

线性回归系数的方差

下面我们来看线性回归系数 β1^\hat{\beta_1}β1^ 和 \hat{\beta_0} 的方差。直观上看，当误差项的方差 σ2\sigma^2σ2 越大时，Var(β1^)\displaystyle \mathrm{Var}( \hat{\beta_1} )Var(β1^) 与 Var(β0^)\displaystyle \mathrm{Var}( \hat{\beta_0} )Var(β0^) 也应该越大。因为误差项的方差越大，随机变量 yiy_iyi 的取值就会越不确定，从而使得线性拟合的直线不确定性越大。

具体地，我们有
Var(β1^)=σ2Sxx\displaystyle \mathrm{Var}( \hat{\beta_1} ) = \frac{ \sigma^2 }{S_{xx}}Var(β1^)=Sxxσ2,
Var(β0^)=σ2(1n+xˉ2Sxx)\displaystyle \mathrm{Var}( \hat{\beta_0} ) = \sigma^2 \left( \frac{ 1 }{ n } + \frac{ \bar{x}^2 }{S_{xx}} \right)Var(β0^)=σ2(n1+Sxxxˉ2)。

对于 β1^\hat{\beta_1}β1^，我们可以直接计算如下，
因为 β1^=∑i=1n(xi−xˉ)Sxx⋅yi\displaystyle \hat{\beta_1} =\sum_{i = 1}^n \frac{ (x_i - \bar{x} ) }{S_{xx} } \cdot y_iβ1^=i=1∑nSxx(xi−xˉ)⋅yi, 我们记 ci=(xi−xˉ)Sxx\displaystyle c_i = \frac{ (x_i - \bar{x} ) }{S_{xx} }ci=Sxx(xi−xˉ)。于是，

Var(β1^)=∑i=1n(ci2⋅Var(yi))=∑i=1n(xi−xˉ)2Sxx2⋅Var(yi)=σ2Sxx2∑i=1n(xi−xˉ)2=σ2Sxx。\begin{aligned} \displaystyle \mathrm{Var}( \hat{\beta_1} ) &= \sum_{i = 1}^n \left( c_i^2 \cdot \mathrm{Var}( y_i ) \right) = \sum_{i = 1}^n \frac{ (x_i - \bar{x} )^2 }{S_{xx}^2 } \cdot \mathrm{Var}( y_i ) \\ &= \frac{\sigma^2}{S_{xx}^2} \sum_{i = 1}^n (x_i - \bar{x} )^2 = \frac{\sigma^2}{S_{xx}}。 \end{aligned}Var(β1^)=i=1∑n(ci2⋅Var(yi))=i=1∑nSxx2(xi−xˉ)2⋅Var(yi)=Sxx2σ2i=1∑n(xi−xˉ)2=Sxxσ2。

下面我们来看 Var(β0^)\displaystyle \mathrm{Var}( \hat{\beta_0} )Var(β0^)。

因为 β0^=yˉ−xˉ⋅β1^\displaystyle \hat{\beta_0} = \bar{y} - \bar{x} \cdot \hat{\beta_1}β0^=yˉ−xˉ⋅β1^，我们有 Var(β0^)=Var(yˉ)+Var(xˉ⋅β1^)−2Cov(yˉ,xˉ⋅β1^)\displaystyle \mathrm{Var}( \hat{\beta_0} ) = \mathrm{Var} (\bar{y}) + \mathrm{Var} ( \bar{x} \cdot \hat{\beta_1}) - 2 \mathrm{Cov} (\bar{y}, \bar{x} \cdot \hat{\beta_1})Var(β0^)=Var(yˉ)+Var(xˉ⋅β1^)−2Cov(yˉ,xˉ⋅β1^)。

我们知道，Var(yˉ)=1n2∑i=1nVar(yi)=1n2∑i=1nVar(ϵi)=σ2n\displaystyle \mathrm{Var} (\bar{y}) = \frac{1}{n^2} \sum_{i=1}^n \mathrm{Var} (y_i) = \frac{1}{n^2} \sum_{i=1}^n \mathrm{Var} (\epsilon_i) = \frac{\sigma^2}{n}Var(yˉ)=n21i=1∑nVar(yi)=n21i=1∑nVar(ϵi)=nσ2。
而 Var(xˉ⋅β1^)=xˉ2⋅Var(β1^)=xˉ2⋅σ2Sxx\displaystyle \mathrm{Var} (\bar{x} \cdot \hat{\beta_1}) = \bar{x}^2 \cdot \mathrm{Var}(\hat{\beta_1}) = \bar{x}^2 \cdot \frac{\sigma^2}{S_{xx}}Var(xˉ⋅β1^)=xˉ2⋅Var(β1^)=xˉ2⋅Sxxσ2。
下面我们来看协方差项 Cov(yˉ,xˉ⋅β1^)=xˉ⋅Cov(yˉ,β1^)\displaystyle \mathrm{Cov} (\bar{y}, \bar{x} \cdot \hat{\beta_1}) = \bar{x} \cdot \mathrm{Cov}(\bar{y}, \hat{\beta_1})Cov(yˉ,xˉ⋅β1^)=xˉ⋅Cov(yˉ,β1^)。

Cov(yˉ,β1^)=Cov(1n∑i=1nyi,∑i=1nxi−xˉSxxyi)=1n∑i=1nCov(yi,xi−xˉSxxyi)=1n∑i=1nxi−xˉSxx⋅σ2=σ2nxi−xˉSxx=0。\begin{aligned} \displaystyle \mathrm{Cov}(\bar{y}, \hat{\beta_1} ) &= \mathrm{Cov} \Big( \frac{1}{n} \sum_{i=1}^n y_i, \sum_{i=1}^n \frac{x_i - \bar{x}}{S_{xx}} y_i \Big) \\ &= \frac{1}{n} \sum_{i=1}^n \mathrm{Cov} \Big( y_i, \frac{x_i - \bar{x}}{S_{xx}} y_i \Big) \\ &= \frac{1}{n} \sum_{i=1}^n \frac{x_i - \bar{x}}{S_{xx}} \cdot \sigma^2 \\ &= \frac{\sigma^2}{n} \frac{x_i - \bar{x}}{S_{xx}} = 0。 \end{aligned}Cov(yˉ,β1^)=Cov(n1i=1∑nyi,i=1∑nSxxxi−xˉyi)=n1i=1∑nCov(yi,Sxxxi−xˉyi)=n1i=1∑nSxxxi−xˉ⋅σ2=nσ2Sxxxi−xˉ=0。
于是，我们有，

Var(β0^)=Var(yˉ)+Var(xˉ⋅β1^)−2Cov(yˉ,xˉ⋅β1^)=σ2n+xˉ2⋅σ2Sxx=σ2(1n+xˉ2Sxx)\displaystyle \mathrm{Var}( \hat{\beta_0} ) = \mathrm{Var} (\bar{y}) + \mathrm{Var} ( \bar{x} \cdot \hat{\beta_1}) - 2 \mathrm{Cov} (\bar{y}, \bar{x} \cdot \hat{\beta_1}) = \frac{\sigma^2}{n} + \bar{x}^2 \cdot \frac{\sigma^2}{S_{xx}} = \sigma^2 \left( \frac{ 1 }{ n } + \frac{ \bar{x}^2 }{S_{xx}} \right)Var(β0^)=Var(yˉ)+Var(xˉ⋅β1^)−2Cov(yˉ,xˉ⋅β1^)=nσ2+xˉ2⋅Sxxσ2=σ2(n1+Sxxxˉ2)。

其余的几个性质

残差项之和为0

我们记残差项为 eie_iei，ei=yi−yi^\displaystyle e_i = y_i - \hat{y_i}ei=yi−yi^。即 ∑i=1n(yi−yi^)=0\displaystyle \sum_{i = 1}^n (y_i - \hat{y_i}) = 0i=1∑n(yi−yi^)=0。

这里我们可以用求 β1^\hat{\beta_1}β1^ 和 β0^\hat{\beta_0}β0^ 的公式来证明。

在文章单变量线性回归的最小二乘法公式中，我们提到在用偏导数求 β1^\hat{\beta_1}β1^ 和 β0^\hat{\beta_0}β0^ 时，我们有β1^\hat{\beta_1}β1^ 和 β0^\hat{\beta_0}β0^ 的表达式如下：

{∂RSS∂β0=2nβ0+2∑i=1nxiβ1−2∑i=1nyi=0∂RSS∂β1=2∑i=1nxi2β1+2∑i=1nxiβ0−2∑i=1nxiyi=0\displaystyle \begin{cases} &\displaystyle \frac{\partial \text{ RSS}}{\partial \beta_0} = 2 n \beta_0 + 2 \sum_{i = 1}^n x_i \beta_1 - 2 \sum_{i = 1}^n y_i = 0 \\ \\ &\displaystyle \frac{\partial \text{ RSS}}{\partial \beta_1} = 2 \sum_{i = 1}^n x_i^2 \beta_1 + 2 \sum_{i = 1}^n x_i \beta_0 - 2 \sum_{i = 1}^n x_i y_i = 0 \\ \end{cases} ⎩⎨⎧∂β0∂ RSS=2nβ0+2i=1∑nxiβ1−2i=1∑nyi=0∂β1∂ RSS=2i=1∑nxi2β1+2i=1∑nxiβ0−2i=1∑nxiyi=0

根据第一个式子，我们把 β1^\hat{\beta_1}β1^ 和 β0^\hat{\beta_0}β0^ 代入，我们有

nβ0^+∑i=1nxiβ1^−∑i=1nyi=0\displaystyle n \hat{\beta_0} + \sum_{i = 1}^n x_i \hat{\beta_1} - \sum_{i = 1}^n y_i = 0nβ0^+i=1∑nxiβ1^−i=1∑nyi=0。
即，∑i=1n(β0^+xiβ1^)−∑i=1nyi=0\displaystyle \sum_{i = 1}^n \left( \hat{\beta_0} + x_i \hat{\beta_1} \right) - \sum_{i = 1}^n y_i = 0i=1∑n(β0^+xiβ1^)−i=1∑nyi=0，亦 ∑i=1n(yi−yi^)=0\displaystyle \sum_{i = 1}^n (y_i - \hat{y_i}) = 0i=1∑n(yi−yi^)=0。

线性拟合直线总会经过 (xˉ,yˉ)(\bar{x}, \bar{y})(xˉ,yˉ) 这个点

拟合直线为 y=β0^+β1^⋅x\displaystyle y = \hat{\beta_0} + \hat{\beta_1} \cdot xy=β0^+β1^⋅x。而我们有
β0^=yˉ−β1^xˉ\displaystyle \hat{\beta_0} = \bar{y} - \hat{\beta_1} \bar{x}β0^=yˉ−β1^xˉ，所以线性拟合直线总会经过 (xˉ,yˉ)(\bar{x}, \bar{y})(xˉ,yˉ) 这个点。

在 xix_ixi 权重下，残差和为0，即∑i=1nxiei=0\displaystyle \sum_{i = 1}^n x_i e_i = 0i=1∑nxiei=0

我们可以直接把 ei=yi−yi^e_i = y_i - \hat{y_i}ei=yi−yi^ 代入。注意到 yi^=β0^+β1^⋅xi\displaystyle \hat{y_i} = \hat{\beta_0} + \hat{\beta_1} \cdot x_iyi^=β0^+β1^⋅xi，我们有

∑i=1nxiei=∑i=1nxi⋅(yi−yi^)=∑i=1nxi⋅(yi−β0^−β1^⋅xi)=∑i=1nxiyi−β0^∑i=1nxi−β1^∑i=1nxi2=∑i=1nxiyi−β1^∑i=1nxi2−(yˉ−β1^⋅xˉ)⋅∑i=1nxi=∑i=1nxiyi−β1^∑i=1nxi2−nxˉyˉ+nβ1^xˉ2=(∑i=1nxiyi−nxˉyˉ)−(∑i=1nxi2−nxˉ2)β1^=Sxy−Sxxβ1^=0\begin{aligned} \displaystyle \sum_{i = 1}^n x_i e_i &= \sum_{i = 1}^n x_i \cdot (y_i - \hat{y_i}) = \sum_{i = 1}^n x_i \cdot (y_i - \hat{\beta_0} - \hat{\beta_1} \cdot x_i) \\ &= \sum_{i = 1}^n x_i y_i - \hat{\beta_0} \sum_{i = 1}^n x_i - \hat{\beta_1} \sum_{i = 1}^n x_i^2 \\ &= \sum_{i = 1}^n x_i y_i - \hat{\beta_1} \sum_{i = 1}^n x_i^2 - (\bar{y} - \hat{\beta_1} \cdot \bar{x} ) \cdot \sum_{i = 1}^n x_i \\ &= \sum_{i = 1}^n x_i y_i - \hat{\beta_1} \sum_{i = 1}^n x_i^2 - n \bar{x} \bar{y} + n \hat{\beta_1} \bar{x}^2 \\ &= \left( \sum_{i = 1}^n x_i y_i - n \bar{x} \bar{y} \right) - \left(\sum_{i = 1}^n x_i^2 - n \bar{x}^2 \right) \hat{\beta_1} \\ &= S_{xy} - S_{xx} \hat{\beta_1} \\ &= 0 \end{aligned}i=1∑nxiei=i=1∑nxi⋅(yi−yi^)=i=1∑nxi⋅(yi−β0^−β1^⋅xi)=i=1∑nxiyi−β0^i=1∑nxi−β1^i=1∑nxi2=i=1∑nxiyi−β1^i=1∑nxi2−(yˉ−β1^⋅xˉ)⋅i=1∑nxi=i=1∑nxiyi−β1^i=1∑nxi2−nxˉyˉ+nβ1^xˉ2=(i=1∑nxiyi−nxˉyˉ)−(i=1∑nxi2−nxˉ2)β1^=Sxy−Sxxβ1^=0

故得证。

在 yi^\hat{y_i}yi^ 权重下，残差和为0，即∑i=1nyi^ei=0\displaystyle \sum_{i = 1}^n \hat{y_i} e_i = 0i=1∑nyi^ei=0

这里我们只须要利用上一个公式，即 ∑i=1nxiei=0\displaystyle \sum_{i = 1}^n x_i e_i = 0i=1∑nxiei=0 即可。

因为我们有 ∑i=1nyi^ei=∑i=1n(β0^+β1^⋅xi)ei=β0^∑i=1nei+β1^∑i=1nxiei=0\displaystyle \sum_{i = 1}^n \hat{y_i} e_i = \sum_{i = 1}^n \big(\hat{\beta_0} + \hat{\beta_1} \cdot x_i \big) e_i = \hat{\beta_0} \sum_{i = 1}^n e_i + \hat{\beta_1} \sum_{i = 1}^n x_i e_i = 0i=1∑nyi^ei=i=1∑n(β0^+β1^⋅xi)ei=β0^i=1∑nei+β1^i=1∑nxiei=0。

模拟

最后我们用 python 程序来模拟“证明” β1^\hat{\beta_1}β1^ 和 β0^\hat{\beta_0}β0^ 分别时 β1\beta_1β1 和 β0\beta_0β0 的无偏估计。

class unbiased_beta:def __init__(self, arr_x: np.array, beta1: float, beta0: float, epsilon: float):#self.N = Nself.X = arr_xself.beta1 = beta1self.beta0 = beta0self.epsilon = epsilonself.Sxx = ((self.X - self.X.mean()) ** 2).sum()self.X_bar = self.X.mean()def estimate_beta(self, N: int) -> tuple:res_beta1, res_beta0 = [], []for i in range(N):#print(i)cur_error = np.random.normal(0, self.epsilon, arr_x.shape)cur_y = self.beta0 + self.beta1 * self.X + cur_errorcur_y_bar = cur_y.mean()Sxy = ((self.X - self.X.mean()) * (cur_y - cur_y_bar)).sum()cur_beta1 = Sxy / self.Sxxcur_beta0 = cur_y_bar - cur_beta1 * self.X_barres_beta1.append(cur_beta1)res_beta0.append(cur_beta0)return np.mean(res_beta1), np.mean(res_beta0)

arr_x = np.array(range(1, 11))
a = unbiased_beta(arr_x, 2, 3, 1)
res = a.estimate_beta(10 ** 5)

res

(1.9988026861047237, 3.0029188805679303)

可以发现，在经过多次的实验之后，我们得到的 β1^\hat{\beta_1}β1^ 和 β0^\hat{\beta_0}β0^ 的平均值是非常接近真实值 β1\beta_1β1 和 β0\beta_0β0 的。

plt.figure(figsize=(8, 6), dpi=100)
plt.hist(res[0], bins=50, density=True);
line_vert = [[2, c] for c in np.linspace(0, 4, 100)]
plt.plot([c[0] for c in line_vert], [c[1] for c in line_vert], '-', linewidth=4)
plt.xlabel("estimated beta1 value", fontsize=20)
plt.ylabel("count", fontsize=20)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)

多元线性回归的系数也是无偏估计

上文中我们介绍了单变量线性回归中，β1^\hat{\beta_1}β1^ 与 β0^\hat{\beta_0}β0^ 是 β1\beta_1β1 与 β0\beta_0β0 的无偏估计 (unbiased) 。下面我们来看多变量线性回归中系数估计的无偏性质。

根据多变量线性回归的 normal equation，我们有

β^=(XTX)−1XTy\displaystyle \mathbf{\hat{\beta}} = (X^TX)^{-1}X^T \mathbf{y}β^=(XTX)−1XTy。

注意到，y=Xβ+ϵ\displaystyle \mathbf{y}= X \mathbf{\beta} + \mathbf{\epsilon}y=Xβ+ϵ。
代入 β^\displaystyle \mathbf{\hat{\beta}}β^ 的表达式，我们有

E[β^]=E[(XTX)−1XTy]=E[(XTX)−1XT(Xβ+ϵ)]=E[(XTX)−1XTXβ]+E[(XTX)−1XTϵ]=β\begin{aligned} \mathbb{E} [ \mathbf{\hat{\beta}} ] &= \mathbb{E} \big[ (X^TX)^{-1}X^T \mathbf{y} \big] = \mathbb{E} \big[ (X^TX)^{-1}X^T ( X \mathbf{\beta} + \mathbf{\epsilon} ) \big] \\ &= \mathbb{E} \big[ (X^TX)^{-1} X^T X \mathbf{\beta} \big] + \mathbb{E} \big[ (X^TX)^{-1} X^T \mathbf{\epsilon} \big] \\ &= \mathbf{\beta} \end{aligned}E[β^]=E[(XTX)−1XTy]=E[(XTX)−1XT(Xβ+ϵ)]=E[(XTX)−1XTXβ]+E[(XTX)−1XTϵ]=β。

于是，我们知道，对于多元线性回归，β^\displaystyle \mathbf{\hat{\beta}}β^ 也是 β\mathbf{\beta}β 的无偏估计。

我们用程序验证如下。

class unbiased_beta_mv:def __init__(self, arr_x: np.array, beta: np.array, epsilon: float):#self.N = Nself.X = arr_xself.beta = betaself.epsilon = epsilondef estimate_beta(self, N: int) -> tuple:res_beta = []for i in range(N):#print(i)cur_error = np.random.normal(0, self.epsilon, arr_x.shape[0])cur_y = self.X @ self.beta + cur_error.reshape((-1, 1))cur_beta_hat = np.linalg.inv(self.X.T @ self.X) @ self.X.T @ cur_yres_beta.append(cur_beta_hat)return np.hstack(res_beta)

beta = np.array([[1], [2], [3], [4]])
arr_x = np.random.normal(0, 1, (10, 3))
arr_x = np.hstack((np.ones((10, 1)), arr_x))

a = unbiased_beta_mv(arr_x, beta, 1)

res = a.estimate_beta(10 ** 5)

np.mean(res, axis=1)

array([1.00255519, 1.99960968, 2.99919333, 4.00008326])

参考文献

[1] Introduction to linear regression analysis, Douglas C. Montgomery, Elizabeth A. Peck, G. Geoffrey Vining, John Wiley & Sons (2021)