在很多研究问题中，我们都希望建立一个定量的解释框架，用以描述一个（或一组）自变量与因变量之间的关系。一个直观的想法是，能否将这种关系表达为一个模型，类似于物理学中使用牛顿第二定律F = m a F=ma F = ma
针对我们感兴趣的因变量（通常称为响应变量），我们希望构建一个关于自变量（通常称为解释变量或预测变量）的模型。然而，由于现实世界的复杂性，我们往往无法完全确定地描述所有影响因变量的因素，因此，一个精确的函数关系很可能是不存在的 。
为了处理这种不确定性，我们允许模型中存在误差项 ，这些误差项代表了我们未能观测或未能纳入模型的随机因素。通过引入误差项，我们得到的模型就会是一个概率模型 。
为了找到最佳拟合数据的模型，我们需要使用最优化 方法。这通常涉及到最小化误差项的某种度量。
上述的逻辑过程就是回归分析的主要内容。出于简便的考虑，回归分析最常用的模型是线性模型，这就是线性回归（linear regression） 。
在行文步骤上，先使用最易于理解的一元情形作为引入，随后通过加入额外特征的需求来介绍多元情形，最后聚焦于回归系数以及线性回归的注意事项。
高斯-马尔可夫定理及其证明
回归模型的矩阵表示
从OLS导出正规方程
记号约定
Symbol Definition Y \mathrm{Y} Y 因变量。一个n n n y i , i = 1 , … , n y_i,i=1,\dots,n y i , i = 1 , … , n i i i Y ^ \hat{\mathrm{Y}} Y ^ X \mathrm{X} X X \mathrm{X} X n × p n \times p n × p X \mathrm{X} X n × ( p + 1 ) n \times (p+1) n × ( p + 1 ) n n n X \mathrm{X} X p p p β \beta β 一个p p p i i i i = 1 , … , p i=1,\dots,p i = 1 , … , p X \mathrm{X} X n × ( p + 1 ) n\times(p+1) n × ( p + 1 ) β 0 \beta_0 β 0 p p p ε \varepsilon ε 一个n n n i i i i = 1 , … , n i=1,\dots,n i = 1 , … , n ε = Y − Y ^ \varepsilon=\mathrm{Y}-\hat{\mathrm{Y}} ε = Y − Y ^
一元自变量下的情形
如果模型只包含1个自变量，那么此时的线性回归可以用
y = k x + b + ε y=kx+b+\varepsilon y = k x + b + ε 来表示。其中k k k x x x y y y b b b k x kx k x x x x y y y
最小二乘法（Ordinary Least Square）
如何求解k k k b b b
y i = k x i + b + ε i ⇔ ε i = y i − ( k x i + b ) . y_i=kx_i+b+\varepsilon_i\Leftrightarrow\varepsilon_i=y_i-(kx_i+b). y i = k x i + b + ε i ⇔ ε i = y i − ( k x i + b ) . 一个很自然的想法是，既然误差是我们被迫引入的，能不能让它尽可能地小？以至于让它接近0呢？
那就是在ε 2 \varepsilon^2 ε 2 k k k b b b ε i 2 \varepsilon_i^2 ε i 2
min k , b ∑ i n ε i 2 = min k , b ∑ i n ( y i − ( k x i + b ) ) 2 利用上面的等式 = min k , b ∑ i n y i 2 − 2 y i ( k x i + b ) + ( k x i + b ) 2 = min k , b ∑ i n − 2 y i ( k x i + b ) + ( k x i + b ) 2 ∑ i n y i 2 与优化目标无关 \begin{aligned}
\min_{k,b}\sum_i^n\varepsilon_i^2&=\min_{k,b}\sum_{i}^n\left(y_i-(kx_i+b)\right)^2&\text{利用上面的等式}\\
&=\min_{k,b}\sum_{i}^ny_i^2-2y_i(kx_i+b)+(kx_i+b)^2&\\
&=\min_{k,b}\sum_{i}^n-2y_i(kx_i+b)+(kx_i+b)^2&\sum_i^ny_i^2\text{与优化目标无关}\\
\end{aligned} k , b min i ∑ n ε i 2 = k , b min i ∑ n ( y i − ( k x i + b ) ) 2 = k , b min i ∑ n y i 2 − 2 y i ( k x i + b ) + ( k x i + b ) 2 = k , b min i ∑ n − 2 y i ( k x i + b ) + ( k x i + b ) 2 利用上面的等式 i ∑ n y i 2 与优化目标无关 将优化对象记为f f f k k k b b b
{ ∂ f ∂ k = ∑ i n − 2 x i y i + 2 b x i + 2 k x i 2 = 0 ∂ f ∂ b = ∑ i n − 2 y i + 2 k x i + 2 b = 0 \left\{\begin{aligned}
\frac{\partial f}{\partial k}&=\sum_i^n-2x_iy_i+2bx_i+2kx_i^2=0\\
\frac{\partial f}{\partial b}&=\sum_i^n-2y_i+2kx_i+2b=0
\end{aligned}\right. ⎩ ⎨ ⎧ ∂ k ∂ f ∂ b ∂ f = i ∑ n − 2 x i y i + 2 b x i + 2 k x i 2 = 0 = i ∑ n − 2 y i + 2 k x i + 2 b = 0 化简后可以得到：
{ ∑ i n − x i y i + b x i + k x i 2 = 0 ⇔ ∑ i n x i ε i = 0 ∑ i n − y i + k x i + b = 0 ⇔ ∑ i n ε i = 0 \left\{\begin{aligned}
&\sum_i^n-x_iy_i+bx_i+kx_i^2=0\Leftrightarrow\sum_i^nx_i\varepsilon_i=0\\
&\sum_i^n-y_i+kx_i+b=0\Leftrightarrow\sum_i^n\varepsilon_i=0
\end{aligned}\right. ⎩ ⎨ ⎧ i ∑ n − x i y i + b x i + k x i 2 = 0 ⇔ i ∑ n x i ε i = 0 i ∑ n − y i + k x i + b = 0 ⇔ i ∑ n ε i = 0 从中可以发现几个有用的等式：
− ∑ i n x i y i + b ∑ i n x i + k ∑ i n x i 2 = 0 ⇔ x y ‾ = k x 2 ‾ + b x ˉ − ∑ i n y i + k ∑ i n x i + n b = 0 ⇔ y ˉ = k x ˉ + b \begin{aligned}
&-\sum_i^nx_iy_i+b\sum_i^nx_i+k\sum_i^nx_i^2=0\Leftrightarrow\overline{xy}=k\overline{x^2}+b\bar{x}\\
&-\sum_i^ny_i+k\sum_i^nx_i+nb=0\Leftrightarrow\bar{y}=k\bar{x}+b
\end{aligned} − i ∑ n x i y i + b i ∑ n x i + k i ∑ n x i 2 = 0 ⇔ x y = k x 2 + b x ˉ − i ∑ n y i + k i ∑ n x i + nb = 0 ⇔ y ˉ = k x ˉ + b 那么就可以解出：
{ k ^ = x y ‾ − x ˉ y ˉ x 2 ‾ − x ˉ 2 这个表示还有更深刻的含义 b ^ = x 2 ‾ y ˉ − x ˉ x y ‾ x 2 ‾ − x ˉ 2 \left\{\begin{aligned}
&\hat{k}=\frac{\overline{xy}-\bar{x}\bar{y}}{\overline{x^2}-\bar{x}^2}&\text{这个表示还有更深刻的含义}\\
&\hat{b}=\frac{\overline{x^2}\bar{y}-\bar{x}\overline{xy}}{\overline{x^2}-\bar{x}^2}&
\end{aligned}\right. ⎩ ⎨ ⎧ k ^ = x 2 − x ˉ 2 x y − x ˉ y ˉ b ^ = x 2 − x ˉ 2 x 2 y ˉ − x ˉ x y 这个表示还有更深刻的含义 对上述推导稍加整理，可以发现，在求k ^ \hat{k} k ^ x ˉ , y ˉ \bar{x},\bar{y} x ˉ , y ˉ b = y ˉ − k x ˉ b=\bar{y}-k\bar{x} b = y ˉ − k x ˉ b ^ \hat{b} b ^
验证
k ^ , b ^ \hat{k},\hat{b} k ^ , b ^ 的确使优化对象取到全局最小值
这里还需要验证如此得到的解的确使得f f f k ∗ ≠ k ^ , b ∗ ≠ b ^ k^*\neq \hat{k},b^*\neq \hat{b} k ∗ = k ^ , b ∗ = b ^ k ^ , b ^ \hat{k},\hat{b} k ^ , b ^
高斯-马尔可夫定理（Gauss-Markov Theorem）
还有一个问题没有解决，上述的优化为什么要使用ε 2 \varepsilon^2 ε 2 ε \varepsilon ε ∣ ε ∣ |\varepsilon| ∣ ε ∣
而一个更加本质的原因是，如果进一步要求ε i ∼ N ( 0 , σ 2 ) \varepsilon_i\sim N(0,\sigma^2) ε i ∼ N ( 0 , σ 2 )
y i ∼ N ( k x i + b , σ 2 ) , y_i\sim N(kx_i+b,\sigma^2), y i ∼ N ( k x i + b , σ 2 ) , 于是
f ( y i ) = 1 2 π σ exp ( − ( y i − ( k x i + b ) ) 2 2 σ 2 ) . f(y_i)=\frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{(y_i-(kx_i+b))^2}{2\sigma^2}\right). f ( y i ) = 2 π σ 1 exp ( − 2 σ 2 ( y i − ( k x i + b ) ) 2 ) . 如果用极大似然法求解，那么似然函数
L ( k , b ) = ∏ i = 1 n f ( y i ) = ( 2 π σ 2 ) − n 2 exp ( − 1 2 σ 2 ∑ i = 1 n ( y i − ( k x i + b ) ) 2 ) \begin{aligned}
L(k,b)&=\prod_{i=1}^nf(y_i)\\
&=(2\pi\sigma^2)^{-\frac{n}{2}}\exp\left(-\frac{1}{2\sigma^2}\sum_{i=1}^n\left(y_i-(kx_i+b)\right)^2\right)
\end{aligned} L ( k , b ) = i = 1 ∏ n f ( y i ) = ( 2 π σ 2 ) − 2 n exp ( − 2 σ 2 1 i = 1 ∑ n ( y i − ( k x i + b ) ) 2 ) 这样就会发现max k , b L ( k , b ) = min k , b ∑ i n ε i 2 \max_{k,b}L(k,b)=\min_{k,b}\sum_i^n\varepsilon_i^2 max k , b L ( k , b ) = min k , b ∑ i n ε i 2
这个假设实际上来自于高斯-马尔可夫定理 ，该定理表明，当以下条件成立时，最小二乘估计就是最优线性无偏估计（Best Linear Unbiased Estimator, BLUE）：
（误差零均值）E ( ε i ) = 0 \mathrm{E}(\varepsilon_i)=0 E ( ε i ) = 0
（误差同方差）V a r ( ε i ) = σ 2 < + ∞ \mathrm{Var}(\varepsilon_i)=\sigma^2<+\infin Var ( ε i ) = σ 2 < + ∞
（误差不相关）C o v ( ε i , ε j ) = 0 , ∀ i ≠ j \mathrm{Cov}(\varepsilon_i,\varepsilon_j)=0,\forall i\neq j Cov ( ε i , ε j ) = 0 , ∀ i = j
先需要搞清楚几个术语的意思。“最优”是指具有最小的方差；“线性”是指因变量是参数的线性组合；“无偏”是指参数估计值的期望与参数的真值相同。
但是，使用当前的语言来证明定理是很繁琐的，为此我们不妨先看多元的线性回归应该如何表示。
多元自变量下的情形
假设有数量为n n n p p p x n , p x_{n,p} x n , p
y i = ∑ j = 0 p β j x i , j + ε i , y_i=\sum_{j=0}^p\beta_jx_{i,j}+\varepsilon_i, y i = j = 0 ∑ p β j x i , j + ε i , 这里β 0 \beta_0 β 0 x i , 0 ≡ 1 x_{i,0}\equiv1 x i , 0 ≡ 1
如何求解参数呢？或许可以继续优化ε 2 \varepsilon^2 ε 2
或者，也可以试试将n n n
矩阵表示的线性回归
根据我们开头的记号约定，我们可以用矩阵乘法来表示线性回归模型：
Y = X β + ε , \mathrm{Y}=\mathrm{X}\beta+\varepsilon, Y = X β + ε , 其中
Y = [ y 1 y 2 ⋮ y n ] , X = [ 1 x 1 , 1 x 1 , 2 ⋯ x 1 , p 1 x 2 , 1 x 2 , 2 ⋯ x 2 , p ⋮ ⋮ ⋮ ⋱ ⋮ 1 x n , 1 x n , 2 ⋯ x n , p ] , β = [ β 0 β 1 ⋮ β p ] , ε = [ ε 1 ε 2 ⋮ ε n ] , \mathrm{Y}=
\begin{bmatrix}
y_1\\
y_2\\
\vdots\\
y_n
\end{bmatrix},
\mathrm{X}=
\begin{bmatrix}
1&x_{1,1}&x_{1,2}&\cdots&x_{1,p}\\
1&x_{2,1}&x_{2,2}&\cdots&x_{2,p}\\
\vdots&\vdots&\vdots&\ddots&\vdots\\
1&x_{n,1}&x_{n,2}&\cdots&x_{n,p}
\end{bmatrix},
\mathrm{\beta}=
\begin{bmatrix}
\beta_0\\
\beta_1\\
\vdots\\
\beta_p
\end{bmatrix},
\mathrm{\varepsilon}=
\begin{bmatrix}
\varepsilon_1\\
\varepsilon_2\\
\vdots\\
\varepsilon_n
\end{bmatrix}, Y = y 1 y 2 ⋮ y n , X = 1 1 ⋮ 1 x 1 , 1 x 2 , 1 ⋮ x n , 1 x 1 , 2 x 2 , 2 ⋮ x n , 2 ⋯ ⋯ ⋱ ⋯ x 1 , p x 2 , p ⋮ x n , p , β = β 0 β 1 ⋮ β p , ε = ε 1 ε 2 ⋮ ε n , 每一行正好对应着我们原先采用的表示。
矩阵表示下的OLS
在这种表示下，我们再进行一次OLS。很显然
ε = Y − X β , \varepsilon=\mathrm{Y}-\mathrm{X}\beta, ε = Y − X β , 而我们的目标min ∑ i = 1 n ε i 2 \min\sum_{i=1}^n\varepsilon_i^2 min ∑ i = 1 n ε i 2 min ∣ ∣ ε ∣ ∣ \min||\varepsilon|| min ∣∣ ε ∣∣ min ε T ε \min\varepsilon^T\varepsilon min ε T ε
min ε T ε = min ( Y − X β ) T ( Y − X β ) . \min \varepsilon^T\varepsilon = \min (\mathrm{Y} - \mathrm{X}\beta)^T(\mathrm{Y} - \mathrm{X}\beta). min ε T ε = min ( Y − X β ) T ( Y − X β ) . 将上式展开，我们得到
min ε T ε = min ( Y T Y − Y T X β − β T X T Y + β T X T X β ) . \min \varepsilon^T\varepsilon = \min (\mathrm{Y} ^T \mathrm{Y} - \mathrm{Y} ^T \mathrm{X} \beta - \beta^T \mathrm{X}^T \mathrm{Y} + \beta^T \mathrm{X}^T \mathrm{X} \beta). min ε T ε = min ( Y T Y − Y T X β − β T X T Y + β T X T X β ) . 对上式关于β \beta β
∂ ∂ β ( ε T ε ) = − 2 X T Y + 2 X T X β . \frac{\partial}{\partial \beta} (\varepsilon^T \varepsilon) = -2\mathrm{X}^T \mathrm{Y} + 2\mathrm{X}^T \mathrm{X} \beta. ∂ β ∂ ( ε T ε ) = − 2 X T Y + 2 X T X β . 首先，Y T Y \mathrm{Y} ^T \mathrm{Y} Y T Y β \beta β
第二项
∂ Y T X β ∂ β = X T Y \frac{\partial\mathrm{Y} ^T \mathrm{X} \beta}{\partial\beta}=\mathrm{X}^T\mathrm{Y} ∂ β ∂ Y T X β = X T Y 这是因为Y T X β \mathrm{Y} ^T \mathrm{X} \beta Y T X β 1 × 1 1\times1 1 × 1
∑ i = 1 n ∑ j = 1 p y i x i j β j . \sum_{i=1}^n\sum_{j=1}^py_ix_{ij}\beta_{j}. i = 1 ∑ n j = 1 ∑ p y i x ij β j . β \beta β p p p β \beta β
( ∑ i = 1 n y i x i 1 , ⋯ , ∑ i = 1 n y i x i p ) T , \left(\sum_{i=1}^ny_ix_{i1},\cdots,\sum_{i=1}^ny_ix_{ip}\right)^T, ( i = 1 ∑ n y i x i 1 , ⋯ , i = 1 ∑ n y i x i p ) T , 它恰好是X T Y \mathrm{X}^T\mathrm{Y} X T Y
第三项
∂ β T X T Y ∂ β = X T Y \frac{\partial\beta^T \mathrm{X}^T\mathrm{Y} }{\partial\beta}=\mathrm{X}^T\mathrm{Y} ∂ β ∂ β T X T Y = X T Y 与上面类似，β T X T Y \beta^T \mathrm{X}^T\mathrm{Y} β T X T Y
∑ i = 1 n ∑ j = 1 p y i x i j β j . \sum_{i=1}^n\sum_{j=1}^py_ix_{ij}\beta_{j}. i = 1 ∑ n j = 1 ∑ p y i x ij β j . 那么求导结果就与上一项一致。
最后一项
∂ β T X T X β ∂ β = 2 X T X β \frac{\partial\beta^T \mathrm{X}^T \mathrm{X} \beta}{\partial\beta}=2\mathrm{X}^T\mathrm{X}\beta ∂ β ∂ β T X T X β = 2 X T X β 同样地，β T X T X β \beta^T \mathrm{X}^T \mathrm{X} \beta β T X T X β
∑ i = 1 n ( ∑ j = 1 p x i j β j ) 2 . \sum_{i=1}^n\left(\sum_{j=1}^px_{ij}\beta_j\right)^2. i = 1 ∑ n ( j = 1 ∑ p x ij β j ) 2 . 这时也是逐个对β \beta β
( 2 ∑ i = 1 n ∑ j = 1 p x i 1 x i j β j , ⋯ , 2 ∑ i = 1 n ∑ j = 1 p x i p x i j β j ) T , \left(2\sum_{i=1}^n\sum_{j=1}^px_{i1}x_{ij}\beta_j,\cdots,2\sum_{i=1}^n\sum_{j=1}^px_{ip}x_{ij}\beta_j\right)^T, ( 2 i = 1 ∑ n j = 1 ∑ p x i 1 x ij β j , ⋯ , 2 i = 1 ∑ n j = 1 ∑ p x i p x ij β j ) T , 这恰好是2 X T X β 2\mathrm{X}^T\mathrm{X}\beta 2 X T X β
总结一下，遇到矩阵求导，首先将其写成求和形式，然后逐项求导。
令导数为零，解得
− 2 X T Y + 2 X T X β = 0. -2\mathrm{X}^T \mathrm{Y} + 2\mathrm{X}^T \mathrm{X} \beta = 0. − 2 X T Y + 2 X T X β = 0. 进一步化简，得到最小二乘法的正规方程
X T X β = X T Y . \mathrm{X}^T \mathrm{X} \beta = \mathrm{X}^T \mathrm{Y} . X T X β = X T Y . 最后，我们可以解出β \beta β
β ^ = ( X T X ) − 1 X T Y . \hat{\beta} = (\mathrm{X}^T \mathrm{X})^{-1} \mathrm{X}^T \mathrm{Y} . β ^ = ( X T X ) − 1 X T Y . 证明高斯马尔可夫定理
线性
线性是指参数β \beta β Y \mathrm{Y} Y
无偏性
对β ^ \hat{\beta} β ^
E ( β ^ ) = ( X T X ) − 1 X T E ( Y ) = ( X T X ) − 1 X T E ( X β + ε ) = β + ( X T X ) − 1 X T E ( ε ) = β \begin{aligned}
\mathrm{E}(\hat{\beta})&=(\mathrm{X}^T \mathrm{X})^{-1} \mathrm{X}^T \mathrm{E}(\mathrm{Y})\\
&=(\mathrm{X}^T \mathrm{X})^{-1} \mathrm{X}^T \mathrm{E}(\mathrm{X}\beta+\varepsilon)\\
&=\beta+(\mathrm{X}^T \mathrm{X})^{-1} \mathrm{X}^T \mathrm{E}(\varepsilon)\\
&=\beta
\end{aligned} E ( β ^ ) = ( X T X ) − 1 X T E ( Y ) = ( X T X ) − 1 X T E ( X β + ε ) = β + ( X T X ) − 1 X T E ( ε ) = β 其中用到了误差零均值的假设。
最佳
首先需要写出V a r ( β ^ ) \mathrm{Var}(\hat{\beta}) Var ( β ^ )
V a r ( β ^ ) = E ( ( β ^ − β ) ( β ^ − β ) T ) = E ( ( X T X ) − 1 X T ε ε T X ( X T X ) − 1 ) \begin{aligned}
\mathrm{Var}(\hat{\beta})&=\mathrm{E}((\hat{\beta}-\beta)(\hat{\beta}-\beta)^T)\\
&=\mathrm{E}\left((\mathrm{X}^T \mathrm{X})^{-1} \mathrm{X}^T \varepsilon\varepsilon^T\mathrm{X}(\mathrm{X}^T\mathrm{X})^{-1}\right)
\end{aligned} Var ( β ^ ) = E (( β ^ − β ) ( β ^ − β ) T ) = E ( ( X T X ) − 1 X T ε ε T X ( X T X ) − 1 ) 利用同方差和误差两两不相关的假设，可以将ε ε T \varepsilon\varepsilon^T ε ε T σ 2 I \sigma^2\mathrm{I} σ 2 I
所以
V a r ( β ^ ) = σ 2 ( X T X ) − 1 . \mathrm{Var}(\hat{\beta})=\sigma^2(\mathrm{X}^T\mathrm{X})^{-1}. Var ( β ^ ) = σ 2 ( X T X ) − 1 . 如果存在一个比OLS更优的线性无偏估计，那么其形式为β ~ = C Y \tilde{\beta}=\mathrm{C}\mathrm{Y} β ~ = CY
C = ( X T X ) − 1 X T + D , \mathrm{C}=(\mathrm{X}^T\mathrm{X})^{-1}\mathrm{X}^T+\mathrm{D}, C = ( X T X ) − 1 X T + D , D \mathrm{D} D
因为无偏，所以E ( β ~ ) = β \mathrm{E}(\tilde{\beta})=\beta E ( β ~ ) = β
E ( β ~ ) = E ( C Y ) = E ( ( ( X T X ) − 1 X T + D ) ( X β + ε ) ) = ( I + D X ) β \begin{aligned}
\mathrm{E}(\tilde{\beta})&=\mathrm{E}(\mathrm{C}\mathrm{Y})\\
&=\mathrm{E}\left(((\mathrm{X}^T\mathrm{X})^{-1}\mathrm{X}^T+\mathrm{D})(\mathrm{X}\beta+\varepsilon)\right)\\
&=(\mathrm{I}+\mathrm{D}\mathrm{X})\beta
\end{aligned} E ( β ~ ) = E ( CY ) = E ( (( X T X ) − 1 X T + D ) ( X β + ε ) ) = ( I + DX ) β 知道D X = 0 \mathrm{D}\mathrm{X}=0 DX = 0
而方差
V a r ( β ~ ) = V a r ( C Y ) = σ 2 C C T = σ 2 ( X T X ) − 1 + σ 2 D D T = V a r ( β ^ ) + σ 2 D D T \begin{aligned}
\mathrm{Var}(\tilde{\beta})&=\mathrm{Var}(\mathrm{C}\mathrm{Y})\\
&=\sigma^2\mathrm{C}\mathrm{C}^T\\
&=\sigma^2(\mathrm{X}^T\mathrm{X})^{-1}+\sigma^2\mathrm{D}\mathrm{D}^T\\
&=\mathrm{Var}(\hat{\beta})+\sigma^2\mathrm{D}\mathrm{D}^T
\end{aligned} Var ( β ~ ) = Var ( CY ) = σ 2 C C T = σ 2 ( X T X ) − 1 + σ 2 D D T = Var ( β ^ ) + σ 2 D D T 就不可能小于V a r ( β ^ ) \mathrm{Var}(\hat{\beta}) Var ( β ^ )
关于线性回归的注意事项
本章对线性回归的原理做了基本说明，但是仍然需要强调一些事情：