前面写了leverage 杠杆的计算以及其意义
主要是为后面的内容做一些铺垫.Cook’s distance起源于提出这个名词的统计学家Cook,用于删除一个样本后,对模型的影响。
假设有如下模型
y
=
X
β
+
ϵ
,
X
∈
R
m
×
p
{\mathbf{y}}= {\mathbf{X}}{{\beta}}+\epsilon,X \in \mathbb{R}^{m \times p}
y=Xβ+ϵ,X∈Rm×p
β
^
=
(
X
T
X
)
−
1
X
T
y
⇒
y
^
=
X
β
^
\hat{\beta}= (X^TX)^{-1}X^Ty\Rightarrow \hat{y}=X\hat{\beta}
β^=(XTX)−1XTy⇒y^=Xβ^
X
(
−
i
)
,
y
(
−
i
)
X_{(-i)},y_{(-i)}
X(−i),y(−i)表示从原来数据中去掉第i行数据
β
(
−
i
)
^
=
(
X
(
−
i
)
T
X
(
−
i
)
)
−
1
X
(
−
i
)
T
y
⇒
y
^
(
−
i
)
=
X
β
^
(
−
i
)
\hat{\beta_{(-i)}}=(X_{(-i)}^TX_{(-i)})^{-1}X_{(-i)}^Ty\Rightarrow\hat{y}_{(-i)}=X\hat{\beta}_{(-i)}
β(−i)^=(X(−i)TX(−i))−1X(−i)Ty⇒y^(−i)=Xβ^(−i)
e
=
y
−
y
^
⇒
s
2
=
(
y
−
y
^
)
T
(
y
−
y
^
)
/
(
n
−
p
)
=
e
T
e
n
−
p
e=y-\hat{y}\Rightarrow s^2=(y-\hat{y})^T(y-\hat{y})/(n-p)=\frac{e^Te}{n-p}
e=y−y^⇒s2=(y−y^)T(y−y^)/(n−p)=n−peTe
n-p表示自由度,显然,这个公式不适合n<=p的情况,对于高维的情况可以参考相应的扩展版。
对第i个样本的cook距离表示如下
D
i
=
(
y
^
(
−
i
)
−
y
^
)
T
(
y
^
(
−
i
)
−
y
^
)
p
s
2
=
(
β
(
−
i
)
^
−
β
^
)
T
X
T
X
(
β
(
−
i
)
^
−
β
^
)
p
s
2
D_i=\frac{(\hat{y}_{(-i)}-\hat{y})^T(\hat{y}_{(-i)}-\hat{y})}{ps^2}=\frac{(\hat{\beta_{(-i)}}-\hat{\beta})^TX^TX(\hat{\beta_{(-i)}}-\hat{\beta})}{ps^2}
Di=ps2(y^(−i)−y^)T(y^(−i)−y^)=ps2(β(−i)^−β^)TXTX(β(−i)^−β^)
上式的变量的平方和,让人很容易想起卡方分布
X
2
\mathcal{X^2}
X2。
两个卡方的相除又让人想到方差齐性检测
F
(
p
,
m
−
p
,
1
−
α
)
F(p,m-p,1-\alpha)
F(p,m−p,1−α)分布,这是
D
i
D_i
Di的主要意义所在。利用了分布的概率
D
i
<
=
F
(
p
,
m
−
p
,
1
−
α
)
D_i<=F(p,m-p,1-\alpha)
Di<=F(p,m−p,1−α)去估计样本的异常情况,显然更加科学,有技术含量。
从表面上看,如果要实现这个功能,需要借助留一法去处理,显然这样做会带来很大的运算量,使得算法的实现变得困难。借助以下公式,使得运算简单
β
^
−
β
^
−
i
=
(
X
T
X
)
−
1
x
i
1
−
v
i
(
y
i
−
x
i
T
β
^
)
\hat{\beta}-\hat{\beta}_{{-i}}=\frac{(X^TX)^{-1}x_i}{1-v_i}(y_i-x_i^T\hat{\beta})
β^−β^−i=1−vi(XTX)−1xi(yi−xiTβ^)
这里,
x
i
x_i
xi表示第i个样本,即X的第i行。
v
i
=
x
i
T
(
X
T
X
)
−
1
x
i
v_i=x_i^T(X^TX)^{-1}x_i
vi=xiT(XTX)−1xi
简略证明如下:
我们对X做行交换,y做相应的变换,是不会影响
β
\beta
β的估计。因此,有
X
=
[
X
(
−
i
)
x
i
T
]
,
y
=
[
y
(
−
i
)
y
i
]
X=\begin{bmatrix} X_{(-i)}\\ x_i^T \end{bmatrix},y=\begin{bmatrix} y_{(-i)}\\ y_i \end{bmatrix}
X=[X(−i)xiT],y=[y(−i)yi]
由于
X
=
[
x
1
T
⋯
x
m
T
]
X = \begin{bmatrix} x_1^T\\ \cdots\\ x_m^T \end{bmatrix}
X=⎣⎡x1T⋯xmT⎦⎤,得到
X
T
X
=
[
x
1
,
⋯
,
x
m
]
[
x
1
T
⋯
x
m
T
]
=
∑
i
=
1
m
x
i
x
i
T
=
X
(
−
i
)
T
X
(
−
i
)
+
x
i
x
i
T
X^TX=[x_1,\cdots,x_m]\begin{bmatrix} x_1^T\\ \cdots\\ x_m^T \end{bmatrix}=\sum_{i=1}^mx_ix_i^T=X_{(-i)}^TX_{(-i)}+x_ix_i^T
XTX=[x1,⋯,xm]⎣⎡x1T⋯xmT⎦⎤=i=1∑mxixiT=X(−i)TX(−i)+xixiT
由于
(
A
+
U
V
′
)
−
1
=
A
−
1
−
(
A
−
1
U
V
′
A
−
1
)
/
(
1
+
V
′
A
−
1
U
)
(A + UV')^{-1} = A^{-1} - (A^{-1}UV'A^{-1})/(1 + V'A^{-1}U)
(A+UV′)−1=A−1−(A−1UV′A−1)/(1+V′A−1U)
令
A
=
X
(
−
i
)
T
X
(
−
i
)
A =X_{(-i)}^TX_{(-i)}
A=X(−i)TX(−i)
(
X
T
X
)
−
1
=
(
X
(
−
i
)
T
X
(
−
i
)
+
x
i
x
i
T
)
−
1
=
A
−
1
−
A
−
1
x
i
x
i
T
A
−
1
/
(
1
+
x
i
T
A
−
1
x
i
)
(X^TX)^{-1}=(X_{(-i)}^TX_{(-i)}+x_ix_i^T)^{-1}=A^{-1}-A^{-1}x_ix_i^TA^{-1}/(1+x_i^TA^{-1}x_i)
(XTX)−1=(X(−i)TX(−i)+xixiT)−1=A−1−A−1xixiTA−1/(1+xiTA−1xi)
X
T
y
=
[
X
(
−
i
)
x
i
T
]
T
[
y
(
−
i
)
y
i
]
=
X
(
−
i
)
T
y
(
−
i
)
+
x
i
y
i
X^Ty=\begin{bmatrix} X_{(-i)}\\ x_i^T \end{bmatrix}^T\begin{bmatrix} y_{(-i)}\\ y_i \end{bmatrix}=X_{(-i)}^Ty_{(-i)}+x_iy_i
XTy=[X(−i)xiT]T[y(−i)yi]=X(−i)Ty(−i)+xiyi
令
w
i
=
x
i
T
(
A
)
−
1
x
i
w_{i}=x_i^T(A)^{-1}x_i
wi=xiT(A)−1xi
β
^
=
(
X
T
X
)
−
1
X
T
y
=
A
−
1
X
(
−
i
)
T
y
(
−
i
)
−
A
−
1
x
i
x
i
T
A
−
1
X
(
−
i
)
T
y
(
−
i
)
/
(
1
+
x
i
T
A
−
1
x
i
)
+
A
−
1
x
i
y
i
−
A
−
1
x
i
x
i
T
A
−
1
x
i
y
i
/
(
1
+
x
i
T
A
−
1
x
i
)
=
(
I
−
A
−
1
x
i
x
i
T
/
(
1
+
w
i
)
)
β
(
−
i
)
+
A
−
1
x
i
y
i
/
(
1
+
w
i
)
\hat{\beta}=(X^TX)^{-1}X^Ty=A^{-1}X_{(-i)}^Ty_{(-i)}-A^{-1}x_ix_i^TA^{-1}X_{(-i)}^Ty_{(-i)}/(1+x_i^TA^{-1}x_i)+\\ A^{-1}x_iy_i-A^{-1}x_ix_i^TA^{-1}x_iy_i/(1+x_i^TA^{-1}x_i)\\ =(I-A^{-1}x_ix_i^T/(1+w_i))\beta_{(-i)}+A^{-1}x_iy_i/(1+w_i)\\
β^=(XTX)−1XTy=A−1X(−i)Ty(−i)−A−1xixiTA−1X(−i)Ty(−i)/(1+xiTA−1xi)+A−1xiyi−A−1xixiTA−1xiyi/(1+xiTA−1xi)=(I−A−1xixiT/(1+wi))β(−i)+A−1xiyi/(1+wi)
由此推得
x
i
T
β
^
=
(
x
i
T
−
w
i
x
i
T
/
(
1
+
w
i
)
)
β
^
(
−
i
)
+
w
i
y
i
/
(
1
+
w
i
)
⇒
x
i
T
β
^
=
x
i
T
β
^
(
−
i
)
/
(
1
+
w
i
)
+
y
i
−
y
i
/
(
1
+
w
i
)
⇒
x
i
T
β
^
−
y
i
=
(
x
i
T
β
^
(
−
i
)
−
y
i
)
/
(
1
+
w
i
)
x_i^T\hat{\beta}=(x_i^T-w_ix_i^T/(1+w_i))\hat{\beta}_{(-i)}+w_iy_i/(1+w_i)\Rightarrow\\ x_i^T\hat{\beta}=x_i^T\hat{\beta}_{(-i)}/(1+w_i)+y_i-y_i/(1+w_i)\Rightarrow\\ x_i^T\hat{\beta}-y_i=(x_i^T\hat{\beta}_{(-i)}-y_i)/(1+w_i)
xiTβ^=(xiT−wixiT/(1+wi))β^(−i)+wiyi/(1+wi)⇒xiTβ^=xiTβ^(−i)/(1+wi)+yi−yi/(1+wi)⇒xiTβ^−yi=(xiTβ^(−i)−yi)/(1+wi)
β ^ − β ^ ( − i ) = A − 1 x i ( y i − x i T β ^ ( − i ) ) / ( 1 + w i ) = A − 1 x i ( y i − x i T β ^ ) \hat{\beta}-\hat{\beta}_{(-i)}=A^{-1}x_i(y_i-x_i^T\hat{\beta}_{(-i)})/(1+w_i)=A^{-1}x_i(y_i-x_i^T\hat{\beta}) β^−β^(−i)=A−1xi(yi−xiTβ^(−i))/(1+wi)=A−1xi(yi−xiTβ^)
由于
X
(
−
i
)
T
X
(
−
i
)
=
X
T
X
−
x
i
x
i
T
X_{(-i)}^TX_{(-i)}=X^TX-x_ix_i^T
X(−i)TX(−i)=XTX−xixiT
(
X
(
−
i
)
T
X
(
−
i
)
)
−
1
=
(
X
T
X
)
−
1
+
(
X
T
X
)
−
1
x
i
x
i
T
(
X
T
X
)
−
1
/
(
1
−
v
i
)
(X_{(-i)}^TX_{(-i)})^{-1}=(X^TX)^{-1}+(X^TX)^{-1}x_ix_i^T(X^TX)^{-1}/(1-v_i)
(X(−i)TX(−i))−1=(XTX)−1+(XTX)−1xixiT(XTX)−1/(1−vi)
v
i
=
x
i
T
X
T
X
x
i
v_i=x_i^TX^TXx_i
vi=xiTXTXxi,可以推得
(
X
(
−
i
)
T
X
(
−
i
)
)
−
1
x
i
=
(
X
T
X
)
−
1
x
i
+
(
X
T
X
)
−
1
x
i
x
i
T
(
X
T
X
)
−
1
x
i
/
(
1
−
v
i
)
=
(
X
T
X
)
−
1
x
i
/
(
1
−
v
i
)
(X_{(-i)}^TX_{(-i)})^{-1}x_i=(X^TX)^{-1}x_i+(X^TX)^{-1}x_ix_i^T(X^TX)^{-1}x_i/(1-v_i)\\ =(X^TX)^{-1}x_i/(1-v_i)
(X(−i)TX(−i))−1xi=(XTX)−1xi+(XTX)−1xixiT(XTX)−1xi/(1−vi)=(XTX)−1xi/(1−vi)
得到
β
^
−
β
^
(
−
i
)
=
(
X
T
X
)
−
1
x
i
1
−
v
i
(
y
i
−
x
i
T
β
^
)
\hat{\beta}-\hat{\beta}_{(-i)}=\frac{(X^TX)^{-1}x_i}{1-v_i}(y_i-x_i^T\hat{\beta})
β^−β^(−i)=1−vi(XTX)−1xi(yi−xiTβ^)
代入
D
i
D_i
Di公式得到
D
i
=
(
y
i
−
x
i
T
β
^
s
1
−
v
i
)
2
v
i
p
(
1
−
v
i
)
D_i = (\frac{y_i-x_i^T\hat{\beta}}{s\sqrt{1-v_i}})^2\frac{v_i}{p(1-v_i)}
Di=(s1−viyi−xiTβ^)2p(1−vi)vi
可以看出
D
i
D_i
Di考虑了样本i的两部分信息,前者是学生化后的残差,后者反应了该样本的杠杆值