异常检测 cook distance

前面写了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β+ϵ,XRm×p
β ^ = ( X T X ) − 1 X T y ⇒ y ^ = X β ^ \hat{\beta}= (X^TX)^{-1}X^Ty\Rightarrow \hat{y}=X\hat{\beta} β^=(XTX)1XTyy^=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)Tyy^(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=yy^s2=(yy^)T(yy^)/(np)=npeTe
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,mp,1α)分布,这是 D i D_i Di的主要意义所在。利用了分布的概率 D i < = F ( p , m − p , 1 − α ) D_i<=F(p,m-p,1-\alpha) Di<=F(p,mp,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=1vi(XTX)1xi(yixiTβ^)
这里, 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=x1TxmT,得到
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]x1TxmT=i=1mxixiT=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=A1(A1UVA1)/(1+VA1U)
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=A1A1xixiTA1/(1+xiTA1xi)

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=A1X(i)Ty(i)A1xixiTA1X(i)Ty(i)/(1+xiTA1xi)+A1xiyiA1xixiTA1xiyi/(1+xiTA1xi)=(IA1xixiT/(1+wi))β(i)+A1xiyi/(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β^=(xiTwixiT/(1+wi))β^(i)+wiyi/(1+wi)xiTβ^=xiTβ^(i)/(1+wi)+yiyi/(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)=A1xi(yixiTβ^(i))/(1+wi)=A1xi(yixiTβ^)

由于 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)=XTXxixiT
( 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/(1vi)
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/(1vi)=(XTX)1xi/(1vi)
得到
β ^ − β ^ ( − 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)=1vi(XTX)1xi(yixiTβ^)
代入 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=(s1vi yixiTβ^)2p(1vi)vi
可以看出 D i D_i Di考虑了样本i的两部分信息,前者是学生化后的残差,后者反应了该样本的杠杆值

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值