最小二乘估计矩阵形式的推导
最近写文章有用到一些算法,自己推一下,顺便mark下来。
这么久我才发现csdn居然都能写Tex了(666)。
考虑一般线性回归模型(OLR)
考虑只含有一个指标的一般线性回归模型(ordinary linear regression model)有如下形式:
y
i
=
β
0
+
β
1
x
i
1
+
ϵ
,
i
=
1
,
2
,
…
,
n
y_i=\beta_0+\beta_1x_{i1}+\epsilon,i=1,2,\dots,n
yi=β0+β1xi1+ϵ,i=1,2,…,n
显然这是基于
n
n
n个观测数据或者叫样本的模型形式。其中
β
0
\beta_0
β0称为截距项系数,
β
1
\beta_1
β1称为
x
1
x_1
x1的回归系数,它们都是未知的常值参数。
ϵ
\epsilon
ϵ是不能被观测到的随机误差项,并且满足
E
(
ϵ
)
=
0
E(\epsilon)=0
E(ϵ)=0,
V
a
r
(
ϵ
)
=
σ
2
>
0
\mathrm{Var(\epsilon)}=\sigma^2>0
Var(ϵ)=σ2>0。其实是有
x
0
x_0
x0的,只是通常认为
x
0
=
1
x_0=1
x0=1。还有一个关键的假设就是
x
x
x不是随机变量(
x
x
x要都随机了,这模型就没法玩了)。
实际上我们所研究的问题往往包含多个指标。那么这些指标
(
x
1
,
x
2
,
.
.
.
,
x
p
)
(x_1,x_2,...,x_p)
(x1,x2,...,xp)就对对应
(
β
0
,
β
1
,
.
.
.
,
β
p
)
(\beta_0,\beta_1,...,\beta_p)
(β0,β1,...,βp)个回归系数,这个时候模型的形式就变成了多元线性回归模型:
y
i
=
β
0
+
β
1
x
i
1
+
β
2
x
i
2
+
⋯
+
β
p
x
i
p
+
ϵ
i
,
i
=
1
,
2
,
…
,
n
y_i=\beta_0+\beta_1x_{i1}+\beta_2x_{i2}+\cdots+\beta_px_{ip}+\epsilon_i , i=1,2,\dots,n
yi=β0+β1xi1+β2xi2+⋯+βpxip+ϵi,i=1,2,…,n
所以为了简化计算和书写方便,我们可以把它写成矩阵的形式:
Y
=
X
β
+
ϵ
Y=X\boldsymbol{\beta}+\boldsymbol{\epsilon}
Y=Xβ+ϵ
Y
=
[
y
1
y
2
⋮
y
n
]
X
=
[
1
x
11
⋯
x
1
p
1
x
21
⋯
x
2
p
⋮
⋮
⋮
⋮
1
x
n
1
⋯
x
n
p
]
β
=
[
β
0
β
1
⋮
β
p
]
ε
=
[
ε
1
ε
2
⋮
ε
n
]
Y=\begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \\ \end{bmatrix} X=\begin{bmatrix} 1 & x_{11} & \cdots & x_{1p} \\ 1 & x_{21} & \cdots & x_{2p} \\ \vdots & \vdots & \vdots &\vdots \\ 1 & x_{n1} & \cdots & x_{np} \\ \end{bmatrix} \boldsymbol{\beta}=\begin{bmatrix} \beta_0\\ \beta_1\\ \vdots\\ \beta_p\\ \end{bmatrix} \boldsymbol{\varepsilon}=\begin{bmatrix} \varepsilon_1\\ \varepsilon_2\\ \vdots\\ \varepsilon_n\\ \end{bmatrix}
Y=⎣⎢⎢⎢⎡y1y2⋮yn⎦⎥⎥⎥⎤X=⎣⎢⎢⎢⎡11⋮1x11x21⋮xn1⋯⋯⋮⋯x1px2p⋮xnp⎦⎥⎥⎥⎤β=⎣⎢⎢⎢⎡β0β1⋮βp⎦⎥⎥⎥⎤ε=⎣⎢⎢⎢⎡ε1ε2⋮εn⎦⎥⎥⎥⎤
其中
X
X
X称为设计矩阵(只是习惯叫法),
Y
Y
Y就不多说了。同样也有一些前提:
X
X
X必须是列满秩;随机误差向量
ε
\boldsymbol{\varepsilon}
ε要满足高斯-马尔科夫条件(1829年,高斯提供了最小二乘法的优化效果强于其他方法的证明,因此被称为高斯-马尔可夫定理):
(i)
E
(
ε
)
=
0
E(\boldsymbol{\varepsilon})=0
E(ε)=0
(ii)
V
a
r
(
ε
)
=
σ
2
I
\mathrm{Var(\boldsymbol{\varepsilon)}}=\sigma^2\boldsymbol{I}
Var(ε)=σ2I
最小二乘估计
最小二乘估计法
(
L
S
E
)
(LSE)
(LSE),它和机器学习领域的梯度下降法还是有一定的区别的(后者没有这么多假设,实用性更广泛),准确的来讲
L
E
S
LES
LES只是一种算法,因为随机误差向量
ϵ
\boldsymbol{\epsilon}
ϵ并不能被观测,所以回归方程不存在解,我们只能尽可能的去接近真实值从而解出全局最优解,即确定一个
β
^
\hat{\boldsymbol{\beta}}
β^使得
ε
=
Y
−
X
β
\boldsymbol{\varepsilon}=Y-X\boldsymbol{\beta}
ε=Y−Xβ各元素的平方和达到最小,可以记为:
Q
(
β
)
=
∑
i
=
1
n
ε
i
2
=
ε
T
ε
=
(
Y
−
X
β
)
T
(
Y
−
X
β
)
=
(
Y
T
Y
−
2
β
T
X
T
Y
+
β
T
X
T
X
β
)
\begin{aligned} Q(\boldsymbol{\beta}) &=\sum_{i=1}^n\varepsilon_i^2\\ &=\boldsymbol{\varepsilon}^T\boldsymbol{\varepsilon}\\ &=(Y-X\boldsymbol{\beta})^T(Y-X\boldsymbol{\beta})\\ &=(Y^TY-2\boldsymbol{\beta}^TX^TY+\boldsymbol{\beta}^TX^TX\boldsymbol{\beta}) \end{aligned}
Q(β)=i=1∑nεi2=εTε=(Y−Xβ)T(Y−Xβ)=(YTY−2βTXTY+βTXTXβ)
令:
∂
Q
(
β
)
∂
β
=
−
2
X
T
Y
+
2
X
T
X
β
=
0
\frac{\partial{Q(\boldsymbol{\beta})}}{\partial\beta}=-2X^TY+2X^TX\boldsymbol{\beta}=0
∂β∂Q(β)=−2XTY+2XTXβ=0
这里需要一些矩阵求导的概念,接下来我们就可以得到一个叫做正规方程 的东西:
X
T
X
β
=
X
T
Y
X^TX\boldsymbol{\beta}=X^TY
XTXβ=XTY
由
r
a
n
k
(
X
T
X
)
=
r
a
n
k
(
X
)
=
p
+
1
\mathrm{rank}(X^TX)=\mathrm{rank}(X)=p+1
rank(XTX)=rank(X)=p+1知
X
T
X
X^TX
XTX是正定矩阵,所以
X
X
X^X
XX存在逆矩阵,那么正规方法就有唯一解了:
β
^
=
(
β
^
0
,
β
^
1
,
⋯
,
β
^
p
)
T
=
(
X
T
X
)
−
1
X
T
Y
\hat{\boldsymbol{\beta}}=(\hat{\beta}_0,\hat{\beta}_1,\cdots,\hat{\beta}_p)^T=(X^TX)^{-1}X^TY
β^=(β^0,β^1,⋯,β^p)T=(XTX)−1XTY
此时
β
\boldsymbol{\beta}
β的估计就得到了,如果再把它带回到模型中去就有:
Y
^
=
X
β
^
=
X
(
X
T
X
)
−
1
X
T
Y
=
S
Y
\hat{Y}=X\hat{\boldsymbol{\beta}}=X(X^TX)^{-1}X^TY=SY
Y^=Xβ^=X(XTX)−1XTY=SY
一般统计学上称
S
S
S是
Y
Y
Y的帽子矩阵,这个称呼是因为有
S
S
S的存在使
Y
Y
Y带上了帽子(总感觉怪怪的?)接下来看残差:
ε
^
=
Y
−
Y
^
=
(
I
−
H
)
Y
\hat{\boldsymbol{\varepsilon}}=Y-\hat{Y}=(I-H)Y
ε^=Y−Y^=(I−H)Y
I
I
I是
n
n
n阶的单位矩阵,显然残差的总和为0,是因为
Q
(
β
)
Q(\boldsymbol{\beta})
Q(β)对截距项求偏导数等于0时:
−
2
∑
i
=
1
n
[
y
i
−
(
β
0
+
∑
i
=
1
p
β
i
x
i
)
]
=
0
-2\sum_{i=1}^n[y_i-(\beta_0+\sum_{i=1}^p\beta_ix_i)]=0
−2i=1∑n[yi−(β0+i=1∑pβixi)]=0
这个式子很明显表达了当存在截距项时,残差和必然为0,这也是为什么200年前拉普拉斯放弃了最小一乘法。也可以证明最小二乘法得到的估计和最大似然估计的结果是相同的,都是无偏估计。关于最小二乘法的BLUE性质不是本文的重点不再赘述。
补充几个推导过程中用到的矩阵求偏导法则
∂
x
T
a
∂
x
=
∂
a
T
x
∂
x
=
a
\frac{\partial x^Ta}{\partial x}=\frac{\partial a^Tx}{\partial x}=a
∂x∂xTa=∂x∂aTx=a
∂
x
T
A
x
∂
x
=
A
x
+
A
T
x
\frac{\partial x^TAx}{\partial x}=Ax+A^Tx
∂x∂xTAx=Ax+ATx
如果
A
A
A是对称的:
A
x
+
A
T
x
=
2
A
x
Ax+A^Tx=2Ax
Ax+ATx=2Ax.
至此推导过程完毕。
参考文献:梅长林,王宁《近代回归分析方法》[M],科学出版社,2012.