Chap 6 Kernel Smoothing Methods
在这一章首先介绍一种回归方法,通过在X不同位置 x 0 x_0 x0的点拟合简单模型 f x 0 f_{x_0} fx0,从而达到在所有取值空间 ℜ p \Re^p ℜp光滑的目的。其中需要使用到smooth的方法,使用kernel,基于两点之间的距离表示 x i x_i xi对 x 0 x_0 x0的影响权重,是一种memory-based的方法。很多时候Kernel还和一个参数 λ \lambda λ有关,用于表示neighborhood的宽度。这样的思路对于核密度估计以及分类都是很有效的。
需要注意的是,这里的核方法核后续介绍SVM等等方法种使用到的是不一样的。
One-Dimensional Kernel Smoothers
回顾另一种memory-based的方法:KNN。KNN实际上可以表示为:
f
^
(
x
)
=
Ave
(
y
i
∣
x
i
∈
N
k
(
x
)
)
\hat{f}(x)=\text{Ave}(y_i|x_i\in N_k(x))
f^(x)=Ave(yi∣xi∈Nk(x))
其中
N
k
(
x
)
N_k(x)
Nk(x)表示x的k近邻集合。易知这种方法得到的函数估计并不是光滑的,因为K近邻集合的变化本身就不是光滑的,这样的不连续性是不必要的,我们往往会假设观测Y是:
Y
=
f
(
X
)
+
ϵ
Y=f(X)+\epsilon
Y=f(X)+ϵ
其中f是光滑的函数,并具有高阶导数。
所以希望可以进行改进,KNN中对K近邻中的点给予相等的权重,K近邻外的点权重为0,现在的改进就是对所有点都给予一个与距离负相关的权重,并且权函数是光滑的,由此得到的拟合函数也会是光滑的。最常见的加权平均形式就是Nadaraya–Watson weighted average:
f
^
(
x
0
)
=
∑
i
=
1
N
K
λ
(
x
0
,
x
i
)
y
i
∑
i
=
1
N
K
λ
(
x
0
,
x
i
)
,
\hat{f}(x_0)=\frac{\sum_{i=1}^NK_\lambda(x_0,x_i)y_i}{\sum_{i=1}^NK_\lambda(x_0,x_i)},
f^(x0)=∑i=1NKλ(x0,xi)∑i=1NKλ(x0,xi)yi,
其中,
K
λ
(
x
0
,
x
)
=
D
(
∣
x
−
x
0
∣
h
λ
(
x
0
)
)
.
K_\lambda(x_0,x)=D\left(\frac{|x-x_0|}{h_\lambda(x_0)}\right).
Kλ(x0,x)=D(hλ(x0)∣x−x0∣).
D(.)是事先给定的与
h
λ
(
.
)
h_\lambda(.)
hλ(.)无关的函数,越靠近0值越小。因此
h
λ
(
.
)
h_\lambda(.)
hλ(.)实际上是对距离的一种scale,常见的选择又两种:
-
h λ ( x ) = λ h_\lambda(x)=\lambda hλ(x)=λ
此时 λ \lambda λ越大,对距离缩减的越多,也就是真正权重大的window越大,这会导致进行加权的点越多,从而方差越小;进行加权的点距离X越远,从而偏差越大。
由于对于真实的函数f进行连续光滑的假设,所以距离 x 0 x_0 x0充分近,可以认为f(.)近似相等,bias越小
一旦 λ \lambda λ取定为某一合适值,bias几乎是不变的,var则需要根据X的分布密度而定(因为此时权重和绝对距离有关,所以权重较大的区域不论 x 0 x_0 x0周围是怎样的,这个区域的大小都是固定的,平均的点都是距离 x 0 x_0 x0比较近的,bias都比较小,而var则和参与主要加权的点的个数有关,所以和落入这个固定大小的区域的点个数有关,和X的分布有关)
-
h k ( x 0 ) = ∣ x 0 − x [ k ] ∣ h_k(x_0)=|x_0-x_{[k]}| hk(x0)=∣x0−x[k]∣
其中 x [ k ] x_{[k]} x[k]是样本点中距离 x 0 x_0 x0第k近的。这种情况下,var是几乎不变的,bias则与X的分布密度有关(因为此时权重和相对距离有关,所以权重较大的区域和 x 0 x_0 x0周围点的分布情况相关,由于是对k近邻的点加权的,所以权重较大的区域可以认为就是K近邻点的集合区域,因此点的个数相对是固定的,方差几乎不变,但是这个区域的大小并不固定,参与主要加权的点距离 x 0 x_0 x0的距离和X的分布有关,因此bias并不是固定的)
-
当X中存在大量相等的情况,使用nearest-neighbors kernel会出现问题,比如由于K近邻只能看到那几个位置的点,没有充分发挥近邻的思想,因此常用的方法是将对应的Y进行平均作为一个点看待,加权的时候权重需要乘以平均的X的个数以体现异质性
常见的核函数 D ( . ) D(.) D(.)有以下几种:
-
Epanechnikov kernel
D ( t ) = { 3 4 ( 1 − t 2 ) if ∣ t ∣ ≤ 1 ; 0 otherwise . D(t)=\left\{\begin{array}{ll}\frac{3}{4}(1-t^2)&\text{if}|t|\leq1;\\0&\text{otherwise}.\end{array}\right. D(t)={43(1−t2)0if∣t∣≤1;otherwise.
具有紧支集,一阶不可导 -
tri-cube function
D ( t ) = { ( 1 − ∣ t ∣ 3 ) 3 if ∣ t ∣ ≤ 1 ; 0 otherwise D(t)=\left\{\begin{array}{ll}(1-|t|^3)^3&\text{if}|t|\leq1;\\0&\text{otherwise}\end{array}\right. D(t)={(1−∣t∣3)30if∣t∣≤1;otherwise
具有紧支集,一阶可导 -
Gaussian kernel
ϕ ( t ) = 1 2 π exp ( − x 2 2 ) \phi(t)=\frac{1}{\sqrt{2\pi}}\exp(-\frac{x^2}{2}) ϕ(t)=2π1exp(−2x2)
不具有紧支集,但是由于3 σ \sigma σ原则,尾部可以近似认为是0,很常用
Local Linear Regression
如果直接使用以上叙述的局部加权进行预测,则预测的bias在边界点较大,因为核函数是对称的;一种减少误差的方法是使用局部线性回归,它可以把误差控制在一阶导数之后。
局部线性回归:
min
α
(
x
0
)
,
β
(
x
0
)
∑
i
=
1
N
K
λ
(
x
0
,
x
i
)
[
y
i
−
α
(
x
0
)
−
β
(
x
0
)
x
i
]
2
.
\begin{aligned}\min_{\alpha(x_0),\beta(x_0)}\sum_{i=1}^NK_\lambda(x_0,x_i)\left[y_i-\alpha(x_0)-\beta(x_0)x_i\right]^2.\end{aligned}
α(x0),β(x0)mini=1∑NKλ(x0,xi)[yi−α(x0)−β(x0)xi]2.
由此对
x
0
x_0
x0处的点进行估计:
f
^
(
x
0
)
=
α
^
(
x
0
)
+
β
^
(
x
0
)
x
0
.
\begin{aligned}\hat{f}(x_0)=\hat{\alpha}(x_0)+\hat{\beta}(x_0)x_0.\end{aligned}
f^(x0)=α^(x0)+β^(x0)x0.
需要注意的是,\textbf{虽然在拟合的时候,对所有的点进行加权最小二乘估计,但是只能用于估计
x
0
x_0
x0的值!}将系数估计值带入:
f
^
(
x
0
)
=
b
(
x
0
)
T
(
B
T
W
(
x
0
)
B
)
−
1
B
T
W
(
x
0
)
y
=
∑
i
=
1
N
l
i
(
x
0
)
y
i
.
\begin{array}{rcl}\hat{f}(x_0)&=&b(x_0)^T(\mathbf{B}^T\mathbf{W}(x_0)\mathbf{B})^{-1}\mathbf{B}^T\mathbf{W}(x_0)\mathbf{y}\\&=&\sum_{i=1}^Nl_i(x_0)y_i.\end{array}
f^(x0)==b(x0)T(BTW(x0)B)−1BTW(x0)y∑i=1Nli(x0)yi.
可以看到还是对
Y
Y
Y的线性加权。其中
l
i
(
x
0
)
l_i(x_0)
li(x0)有时也被称作equivalent kernel.这个权重相对于直接的核函数加权并不一定是对称的,还可以根据数据灵活调整每一个点的权重,因此可以缓解bias大的情况(图中绿色的点表示了equivalent的权重
l
i
(
x
0
)
l_i(x_0)
li(x0))。
Local Polynomial Regression
将局部线性回归继续推广成局部多项式回归(degree=d):
min
α
(
x
0
)
,
β
j
(
x
0
)
,
j
=
1
,
…
,
d
∑
i
=
1
N
K
λ
(
x
0
,
x
i
)
[
y
i
−
α
(
x
0
)
−
∑
j
=
1
d
β
j
(
x
0
)
x
i
j
]
2
\begin{aligned}\min_{\alpha(x_0),\beta_j(x_0),j=1,\dots,d}\sum_{i=1}^NK_\lambda(x_0,x_i)\left[y_i-\alpha(x_0)-\sum_{j=1}^d\beta_j(x_0)x_i^j\right]^2\end{aligned}
α(x0),βj(x0),j=1,…,dmini=1∑NKλ(x0,xi)[yi−α(x0)−j=1∑dβj(x0)xij]2
由此得到
x
0
x_0
x0处的估计值:
f
^
(
x
0
)
=
α
^
(
x
0
)
+
∑
j
=
1
d
β
^
j
(
x
0
)
x
0
j
.
\begin{aligned}\hat{f}(x_0)=\hat{\alpha}(x_0)+\sum_{j=1}^d\hat{\beta}_j(x_0)x_0^j.\end{aligned}
f^(x0)=α^(x0)+j=1∑dβ^j(x0)x0j.,并且我们有:
b
j
(
x
0
)
=
∑
i
=
1
N
(
x
i
−
x
0
)
j
l
i
(
x
0
)
.
b_{j}(x_{0})=\sum_{i=1}^{N}(x_{i}-x_{0})^{j}l_{i}(x_{0}).
bj(x0)=i=1∑N(xi−x0)jli(x0).
b
0
(
x
0
)
=
1
,
b
k
(
x
0
)
=
0
,
k
=
1
,
2
,
.
.
.
,
d
b_0(x_0)=1,b_k(x_0)=0,k=1,2,...,d
b0(x0)=1,bk(x0)=0,k=1,2,...,d
# 证明:
最小化:
∑
i
=
1
N
K
λ
(
x
0
,
x
i
)
(
y
i
−
α
(
x
0
)
−
∑
j
=
1
d
x
j
j
β
j
(
x
0
)
)
2
\sum_{i=1}^NK_\lambda(x_0,x_i)(y_i-\alpha(x_0)-\sum_{j=1}^d x_j^j\beta_j(x_0))^2
i=1∑NKλ(x0,xi)(yi−α(x0)−j=1∑dxjjβj(x0))2
也即加权最小二乘:
(
Y
−
B
β
)
T
W
(
Y
−
B
β
)
(Y-B\beta)^TW(Y-B\beta)
(Y−Bβ)TW(Y−Bβ)
其中
B
=
[
1
,
X
,
X
2
,
.
.
.
,
X
d
]
B=[\mathrm{1},X,X^2,...,X^d]
B=[1,X,X2,...,Xd],
W
=
d
i
a
g
(
K
λ
(
x
0
,
x
i
)
)
W=diag(K_\lambda(x_0,x_i))
W=diag(Kλ(x0,xi)),记
b
(
x
)
=
[
1
,
x
,
x
2
,
.
.
.
,
x
d
]
T
b(x)=[1,x,x^2,...,x^d]^T
b(x)=[1,x,x2,...,xd]T
则易知:
β
^
=
(
B
T
W
B
)
−
1
B
T
W
Y
\hat \beta=(B^TWB)^{-1}B^TWY
β^=(BTWB)−1BTWY
y
^
(
x
0
)
=
β
^
T
b
(
x
0
)
=
b
(
x
0
)
T
(
B
T
W
B
)
−
1
B
T
W
Y
\hat y(x_0)=\hat \beta^T b(x_0)=b(x_0)^T(B^TWB)^{-1}B^TWY
y^(x0)=β^Tb(x0)=b(x0)T(BTWB)−1BTWY
从而
L
(
x
0
)
T
=
b
(
x
0
)
T
(
B
T
W
B
)
−
1
B
T
W
L(x_0)^T=b(x_0)^T(B^TWB)^{-1}B^TW
L(x0)T=b(x0)T(BTWB)−1BTW
而本结论的证明最关键的一点是:
b
(
x
0
)
T
=
b
(
x
0
)
T
(
B
T
W
B
)
−
1
B
T
W
B
=
L
(
x
0
)
T
B
b(x_0)^T=b(x_0)^T(B^TWB)^{-1}B^TWB=L(x_0)^TB
b(x0)T=b(x0)T(BTWB)−1BTWB=L(x0)TB
所以:
b
(
x
0
)
i
=
L
(
x
0
)
T
B
i
b(x_0)_i=L(x_0)^TB_i
b(x0)i=L(x0)TBi
其中
B
i
B_i
Bi表示B的第i个列向量。由此可以得到:
∑
i
=
1
N
x
i
j
L
i
(
x
0
)
=
x
0
j
,
j
=
0
,
1
,
2
,
.
.
.
d
\sum_{i=1}^N x_i^jL_i(x_0)=x_0^j,j=0,1,2,...d
i=1∑NxijLi(x0)=x0j,j=0,1,2,...d
j=0,则
∑
i
=
1
N
L
i
(
x
0
)
=
1
\sum_{i=1}^N L_i(x_0)=1
∑i=1NLi(x0)=1;j>0:
∑
i
=
1
N
(
x
i
−
x
0
)
j
L
i
(
x
0
)
=
∑
i
=
1
N
(
∑
k
=
0
d
(
j
k
)
x
i
k
x
0
j
−
k
(
−
1
)
j
−
k
)
=
∑
k
=
0
j
(
j
k
)
(
−
1
)
j
−
k
x
0
j
−
k
∑
i
=
1
N
x
i
j
L
i
(
x
0
)
=
∑
k
=
0
j
(
j
k
)
(
−
1
)
j
−
k
x
0
j
−
k
x
0
k
=
0
\begin{aligned} \sum_{i=1}^N(x_i-x_0)^jL_i(x_0)=&\sum_{i=1}^N(\sum_{k=0}^d \begin{pmatrix} j\\ k \end{pmatrix}x_i^kx_0^{j-k}(-1)^{j-k})\\ =&\sum_{k=0}^j\begin{pmatrix} j\\ k \end{pmatrix}(-1)^{j-k}x_0^{j-k}\sum_{i=1}^Nx_i^j L_i(x_0)\\ =&\sum_{k=0}^j\begin{pmatrix} j\\ k \end{pmatrix}(-1)^{j-k}x_0^{j-k}x_0^{k}\\ =&0 \end{aligned}
i=1∑N(xi−x0)jLi(x0)====i=1∑N(k=0∑d(jk)xikx0j−k(−1)j−k)k=0∑j(jk)(−1)j−kx0j−ki=1∑NxijLi(x0)k=0∑j(jk)(−1)j−kx0j−kx0k0
现在来分析bias:
E
f
^
(
x
0
)
=
∑
i
=
1
N
l
i
(
x
0
)
f
(
x
i
)
=
f
(
x
0
)
∑
i
=
1
N
l
i
(
x
0
)
+
f
′
(
x
0
)
∑
i
=
1
N
(
x
i
−
x
0
)
l
i
(
x
0
)
+
f
′
′
(
x
0
)
2
∑
i
=
1
N
(
x
i
−
x
0
)
2
l
i
(
x
0
)
+
.
.
.
+
f
(
d
)
(
x
0
)
d
!
∑
i
=
1
N
(
x
i
−
x
0
)
d
l
i
(
x
0
)
+
R
=
f
(
x
0
)
+
R
\begin{aligned} \operatorname{E}\hat{f}(x_{0})& \begin{aligned}=\quad\sum_{i=1}^Nl_i(x_0)f(x_i)\end{aligned} \\ &=\quad f(x_0)\sum_{i=1}^Nl_i(x_0)+f'(x_0)\sum_{i=1}^N(x_i-x_0)l_i(x_0) \\ &+\frac{f^{\prime\prime}(x_{0})}{2}\sum_{i=1}^{N}(x_{i}-x_{0})^{2}l_{i}(x_{0})+...+\frac{f^{(d)}(x_0)}{d!}\sum_{i=1}^N(x_i-x_0)^dl_i(x_0)+R\\ &=f(x_0)+R \end{aligned}
Ef^(x0)=i=1∑Nli(x0)f(xi)=f(x0)i=1∑Nli(x0)+f′(x0)i=1∑N(xi−x0)li(x0)+2f′′(x0)i=1∑N(xi−x0)2li(x0)+...+d!f(d)(x0)i=1∑N(xi−x0)dli(x0)+R=f(x0)+R
其中R是更高阶的余项。根据
f
(
x
)
f(x)
f(x)光滑性的假设,高阶导数的绝对值会越来越小,因此随着d的增加,偏差逐渐减小。但同时,由于:
V
a
r
(
f
^
(
x
0
)
)
=
∑
i
=
1
N
l
i
2
(
x
0
)
σ
2
Var(\hat f(x_0))=\sum_{i=1}^Nl_i^2(x_0)\sigma^2
Var(f^(x0))=i=1∑Nli2(x0)σ2
再由一系列线性代数运算知道
∣
∣
l
(
x
0
)
∣
∣
2
2
||l(x_0)||_2^2
∣∣l(x0)∣∣22随着d的增加而增加,因此随着d的增加,方差越大,这就有一个trade-off。
总的来说,对比局部线性回归和局部多项式回归,我们有一些结论:
-
局部线性回归可以很好地减少边界点处估计的偏差,局部多项式回归虽然也可以进一步减少偏差,但是会增大方差;
-
局部线性回归在曲线的波峰或波谷上拟合效果较差,此时可以考虑使用局部多项式回归
-
局部多项式回归中,更建议使用偶数次的多项式(\emph{这里其实没有太理解为什么,原文说的是"Asymptotic analysis suggest that local polynomials of odd degree dominate those of even degree. This is largely due to the fact that asymptotically the MSE is dominated by boundary effects."})
-
由于更多的时候是关注外推的结果,所以为了减小方差,更常用的是局部线性回归
Selecting with Width of the Kernel
如果选择
h
λ
(
x
)
=
λ
h_\lambda(x)=\lambda
hλ(x)=λ,有一个问题就是选择
λ
\lambda
λ,也就是选择kernel的width,这对于Epanechnikov或tri-cube kernel,影响的是支集的范围,对Gaussian kernel来说,影响的是标准差。对于局部线性回归来说,如果width
→
\to
→ 0,就相当于只能看见当前点,近似于一个线性插值函数;如果width
→
+
∞
\to+\infty
→+∞,就相当于每个点都看到了并且是等权重的,近似于OLSE。类似的,在局部线性回归中:
f
^
(
x
0
)
=
b
(
x
0
)
T
(
B
T
W
(
x
0
)
B
)
−
1
B
T
W
(
x
0
)
y
\hat{f}(x_{0})\quad=\quad b(x_{0})^{T}(\mathbf{B}^{T}\mathbf{W}(x_{0})\mathbf{B})^{-1}\mathbf{B}^{T}\mathbf{W}(x_{0})\mathbf{y}
f^(x0)=b(x0)T(BTW(x0)B)−1BTW(x0)y
其中B是1和X组合的矩阵,
W
(
x
0
)
W(x_0)
W(x0)则考虑到了
λ
\lambda
λ,和定义岭回归自由度一样可以定义局部线性回归自由度为
t
r
a
c
e
(
S
λ
=
∑
i
=
1
N
l
i
(
x
i
)
trace(S_\lambda=\sum_{i=1}^Nl_i(x_i)
trace(Sλ=∑i=1Nli(xi)(事实上,
S
i
j
=
l
i
(
x
j
)
S_{ij}=l_i(x_j)
Sij=li(xj))。
Local Regression in ℜ p \Re^p ℜp
现在推广到X是多元的情况。如果是degree=d的局部多项式回归,那么就需要衍生出可能的交互项和高阶项,这和没有结点的基函数是类似的。再进一步使用最小二乘回归进行系数估计即可。(注:由于计算kernel得到weight的时候需要计算距离,所以需要先标准化到同一量纲)
对于一维的情况都容易出现边界点对整体估计影响大的情况,在高维更是如此(由维度灾难理论,更多的点会集中在边界处)。局部多项式回归调整d可以减少偏差,但是在高于3维的时候,局部线性回归就已经不再适用了。并且也不再适合直接使用点图来进行可视化,更多的是使用condition plots,也就是在控制其余几个变量不同值的时候,针对剩余1个变量的可视化插值图。
Structured Local Regression Models in ℜ p \Re^p ℜp
当维度和样本量的比例并没有很小的时候,local regression并不会有很好的效果(维度灾难),此时可以根据数据做出一些结构化的假设,例如是具有稀疏的特征等等。
Structured Kernels
\textbf{第一种进行结构化调整的方法是直接调整local regression的核函数},一般的核函数对于X的每一维是等权重的,也就是:
K
λ
(
x
0
,
x
)
=
D
(
∣
∣
x
−
x
0
∣
∣
λ
)
,
K_\lambda(x_0,x)=D\left(\frac{||x-x_0||}{\lambda}\right),
Kλ(x0,x)=D(λ∣∣x−x0∣∣),
由于是等权重的所以一般需要先进行标准化。而structured kernels则是直接进行非等权的处理:
K
λ
,
A
(
x
0
,
x
)
=
D
(
(
x
−
x
0
)
T
A
(
x
−
x
0
)
λ
)
.
K_{\lambda,A}(x_{0},x)=D\left(\frac{(x-x_{0})^{T}\mathbf{A}(x-x_{0})}{\lambda}\right).
Kλ,A(x0,x)=D(λ(x−x0)TA(x−x0)).
其中A是一个半正定的矩阵,可以选为
Σ
^
−
1
\hat \Sigma^{-1}
Σ^−1,或者选为对角阵来直接单独调整每一个维度的权重,总之调整A来进行一些结构上的调整是一种方法,但是我们更倾向使用的是下面介绍的方法。
Structured Regression Functions
由于是使用f来拟合
E
(
Y
∣
X
)
=
F
(
X
1
,
X
2
,
.
.
.
,
X
p
)
E(Y|X)=F(X_1,X_2,...,X_p)
E(Y∣X)=F(X1,X2,...,Xp)所以一种很直接的想法是使用ANOVA中主效应、交互效应加和的分解方法:
f
(
X
1
,
X
2
,
…
,
X
p
)
=
α
+
∑
j
g
j
(
X
j
)
+
∑
k
<
ℓ
g
k
ℓ
(
X
k
,
X
ℓ
)
+
⋯
f(X_1,X_2,\ldots,X_p)=\alpha+\sum_jg_j(X_j)+\sum_{k<\ell}g_{k\ell}(X_k,X_\ell)+\cdots
f(X1,X2,…,Xp)=α+j∑gj(Xj)+k<ℓ∑gkℓ(Xk,Xℓ)+⋯
当然其中高阶的交互项并不会纳入考虑。拟合locally版本的
f
(
X
)
=
α
+
∑
j
=
1
p
g
j
(
X
j
)
f(X)=\alpha+\sum_{j=1}^{p}g_{j}(X_{j})
f(X)=α+∑j=1pgj(Xj)一种方法是:假设已知除了第k项之外其余项的估计,则对于残差和第k项进行拟合,再对每一项重复以上过程直到收敛。需要注意的是,每一次拟合的时候都是一个local regression,以上方法对于低维度的ANOVA分解拟合比较常用。
一个structured model很重要的例子是varying coefficient models,假如把整个X分成两部分,
(
X
i
1
,
X
i
2
,
.
.
.
,
X
i
p
)
(X_{i1},X_{i2},...,X_{ip})
(Xi1,Xi2,...,Xip)以及
Z
=
(
X
j
1
,
X
j
2
,
.
.
.
X
j
q
)
Z=(X_{j1},X_{j2},...X_{jq})
Z=(Xj1,Xj2,...Xjq)这两个集合可以有重叠的部分。构建模型:
f
(
X
)
=
α
(
Z
)
+
β
1
(
Z
)
X
i
1
+
⋯
+
β
p
(
Z
)
X
i
p
.
f(X)=\alpha(Z)+\beta_1(Z)X_{i1}+\cdots+\beta_p(Z)X_{ip}.
f(X)=α(Z)+β1(Z)Xi1+⋯+βp(Z)Xip.
例如一元局部线性回归实际上就是Z=X,
f
(
X
)
=
β
0
(
X
)
+
β
1
(
X
)
X
f(X)=\beta_0(X)+\beta_1(X)X
f(X)=β0(X)+β1(X)X。当给定Z的时候,以上模型是一个线性模型,但
β
\beta
β会随着Z变化而变化。这个模型的参数估计一种很自然的想法就是使用局部线性回归:
min
α
(
z
0
)
,
β
(
z
0
)
∑
i
=
1
N
K
λ
(
z
0
,
z
i
)
(
y
i
−
α
(
z
0
)
−
x
1
i
β
1
(
z
0
)
−
⋯
−
x
p
i
β
p
(
z
0
)
)
2
.
\begin{aligned}\min_{\alpha(z_{0}),\beta(z_{0})}\sum_{i=1}^{N}K_{\lambda}(z_{0},z_{i})\left(y_{i}-\alpha(z_{0})-x_{1i}\beta_{1}(z_{0})-\cdots-x_{pi}\beta_{p}(z_{0})\right)^{2}.\end{aligned}
α(z0),β(z0)mini=1∑NKλ(z0,zi)(yi−α(z0)−x1iβ1(z0)−⋯−xpiβp(z0))2.
【注】:从这里就可以看出local方法的优点,并没有直接建模 β . ( Z ) \beta_.(Z) β.(Z)和Z之间的关系,而是使用local的方法,每次拟合的时候given Z从而拟合所有可能的函数
Local Likelihood and Other Models
回顾上一节的方法,也就是使用local regression+varying coefficient,从某种程度上说,任何参数化的模型都可以使用这种方法进行拟合:
- 对数似然函数具有这样的形式:
l ( β ) = ∑ i = 1 N l ( y i , x i T β ) . l(\beta)=\sum_{i=1}^Nl(y_i,x_i^T\beta). l(β)=i=1∑Nl(yi,xiTβ).
其中 θ i = θ ( x i ) = x i T β \theta_i=\theta(x_i)=x_i^T\beta θi=θ(xi)=xiTβ,可以使用local的方法使得模型更加flexible θ ( x 0 ) = x 0 T β ( x 0 ) \theta(x_0)=x_0^T\beta(x_0) θ(x0)=x0Tβ(x0),也就是:
l ( β ( x 0 ) ) = ∑ i = 1 N K λ ( x 0 , x i ) l ( y i , x i T β ( x 0 ) ) . l(\beta(x_0))=\sum_{i=1}^NK_\lambda(x_0,x_i)l(y_i,x_i^T\beta(x_0)). l(β(x0))=i=1∑NKλ(x0,xi)l(yi,xiTβ(x0)).
使用local likelihood相对于类似逻辑回归这种log-linear model来说,将全局的线性缩小到局部的线性,是更加flexible的。
从这个例子也可以看出locally regression+varying coefficient应用的一般形式,首先对于感兴趣的函数使用varying coefficient的线性加和的形式进行建模,再对目标函数局部使用使用核函数(关于coefficient的变量)进行加权最大/小化
-
将上一种方法进行推广:
l ( θ ( z 0 ) ) = ∑ i = 1 N K λ ( z 0 , z i ) l ( y i , η ( x i , θ ( z 0 ) ) ) . l(\theta(z_0))=\sum_{i=1}^NK_\lambda(z_0,z_i)l(y_i,\eta(x_i,\theta{(z_0)})). l(θ(z0))=i=1∑NKλ(z0,zi)l(yi,η(xi,θ(z0))).
也就是没有局限于只能 θ \theta θ是一个线性的形式。 -
自回归模型中 y t = β 0 + β 1 y t − 1 + β 2 y t − 2 + ⋯ + β k y t − k + ε t . y_t=\beta_{0}+\beta_{1}y_{t-1}+\beta_{2}y_{t-2}+\cdots+\beta_{k}y_{t-k}+\varepsilon_{t}. yt=β0+β1yt−1+β2yt−2+⋯+βkyt−k+εt.,也就是 z t = ( y t − 1 , y t − 2 , … , y t − k ) , y t = z t T β + ϵ t z_{t}=(y_{t-1},y_{t-2},\ldots,y_{t-k}),y_t=z_t^T\beta+\epsilon_t zt=(yt−1,yt−2,…,yt−k),yt=ztTβ+ϵt,再使用最小二乘估计进行拟合就行。而如果使用由核函数 K ( z 0 , z t ) K(z_0,z_t) K(z0,zt)加权的最小二乘加权,并辅佐局部最小二乘回归,则可以得到一个根据t变化系数的自回归模型。
-
局部多分类逻辑回归
∑
i
=
1
N
K
λ
(
x
0
,
x
i
)
{
β
g
i
0
(
x
0
)
+
β
g
i
(
x
0
)
T
(
x
i
−
x
0
)
−
log
[
1
+
∑
k
=
1
J
−
1
exp
(
β
k
0
(
x
0
)
+
β
k
(
x
0
)
T
(
x
i
−
x
0
)
)
]
}
\begin{aligned} \begin{aligned}\sum_{i=1}^NK_\lambda(x_0,x_i)\end{aligned}& \begin{cases}\beta_{g_i0}(x_0)+\beta_{g_i}(x_0)^T(x_i-x_0)\end{cases} \\ &\left.-\log\left[1+\sum\limits_{k=1}^{J-1}\exp\left(\beta_{k0}(x_0)+\beta_k(x_0)^T(x_i-x_0)\right)\right]\right\} \end{aligned}
i=1∑NKλ(x0,xi){βgi0(x0)+βgi(x0)T(xi−x0)−log[1+k=1∑J−1exp(βk0(x0)+βk(x0)T(xi−x0))]}
可以看到由于在最小化的损失函数中,已经围绕
x
0
x_0
x0进行了中心化,所以
x
0
x_0
x0处的估计为:
Pr
^
(
G
=
j
∣
X
=
x
0
)
=
e
β
^
j
0
(
x
0
)
1
+
∑
k
=
1
J
−
1
e
β
^
k
0
(
x
0
)
.
\hat{\Pr}(G=j|X=x_0)=\frac{e^{\hat{\beta}_{j0}(x_0)}}{1+\sum_{k=1}^{J-1}e^{\hat{\beta}_{k0}(x_0)}}.
Pr^(G=j∣X=x0)=1+∑k=1J−1eβ^k0(x0)eβ^j0(x0).
这个模型相较于一般的多分类逻辑回归是更加flexible的。
特别地,以上方法和使用Nadaraya-Watson kernel smoother估计得到的概率是完全相等的。 对于Nadaraya-Watson kernel smoother有:
P
(
G
=
j
∣
X
=
x
0
)
=
∑
i
∈
G
j
K
λ
(
x
0
,
x
i
)
∑
i
=
1
N
K
λ
(
x
0
,
x
i
)
∝
∑
i
∈
G
j
K
λ
(
x
0
,
x
i
)
.
\mathrm{P}(G=j|X=x_0)=\frac{\sum_{i\in G_j}K_\lambda(x_0,x_i)}{\sum_{i=1}^NK_\lambda(x_0,x_i)}\propto\sum_{i\in G_j}K_\lambda(x_0,x_i).
P(G=j∣X=x0)=∑i=1NKλ(x0,xi)∑i∈GjKλ(x0,xi)∝i∈Gj∑Kλ(x0,xi).
而对于以上的局部多项式逻辑回归,损失函数对每一个系数求导:
∂
l
(
β
,
x
0
)
∂
β
j
0
=
∑
i
=
1
N
K
λ
(
x
0
,
x
i
)
(
1
(
y
i
=
j
)
−
e
β
j
0
+
β
j
T
(
x
i
−
x
0
)
1
+
∑
l
=
1
K
−
1
e
β
l
0
+
β
l
T
(
x
i
−
x
0
)
)
=
∑
i
∈
G
j
K
λ
(
x
0
,
x
i
)
−
e
β
j
0
⋅
∑
i
=
1
N
K
λ
(
x
0
,
x
i
)
e
β
j
T
(
x
i
−
x
0
)
1
+
∑
l
=
1
K
−
1
e
β
l
0
+
β
l
T
(
x
i
−
x
0
)
=
0.
\begin{aligned} \frac{\partial l(\beta,x_{0})}{\partial\beta_{j0}}& =\sum_{i=1}^NK_\lambda(x_0,x_i)\left(\mathbf{1}(y_i=j)-\frac{e^{\beta_{j0}+\beta_j^T(x_i-x_0)}}{1+\sum_{l=1}^{K-1}e^{\beta_{l0}+\beta_l^T(x_i-x_0)}}\right) \\ &=\sum_{i\in G_j}K_\lambda(x_0,x_i)-e^{\beta_{j0}}\cdot\sum_{i=1}^NK_\lambda(x_0,x_i)\frac{e^{\beta_j^T(x_i-x_0)}}{1+\sum_{l=1}^{K-1}e^{\beta_{l0}+\beta_l^T(x_i-x_0)}} \\ &=0. \end{aligned}
∂βj0∂l(β,x0)=i=1∑NKλ(x0,xi)(1(yi=j)−1+∑l=1K−1eβl0+βlT(xi−x0)eβj0+βjT(xi−x0))=i∈Gj∑Kλ(x0,xi)−eβj0⋅i=1∑NKλ(x0,xi)1+∑l=1K−1eβl0+βlT(xi−x0)eβjT(xi−x0)=0.
从而有:
Pr
^
(
G
=
j
∣
X
=
x
0
)
∝
exp
(
β
^
j
0
)
∝
∑
i
∈
G
j
K
λ
(
x
0
,
x
i
)
.
\hat{\Pr}(G=j|X=x_0)\propto\exp(\hat{\beta}_{j0})\propto\sum_{i\in G_j}K_\lambda(x_0,x_i).
Pr^(G=j∣X=x0)∝exp(β^j0)∝i∈Gj∑Kλ(x0,xi).
也就是二者是等价的!
Kernel Density Estimation and Classification
本节介绍核密度估计以及由此引申出的分类方法。
Kernel Density Estimation
一种最容易使用的密度估计方法(直方图法)是:
f
^
X
(
x
0
)
=
#
x
i
∈
N
(
x
0
)
N
λ
,
\hat{f}_X(x_0)=\frac{\#x_i\in\mathcal{N}(x_0)}{N\lambda},
f^X(x0)=Nλ#xi∈N(x0),
其中
N
(
x
0
)
N(x_0)
N(x0)表示以
x
0
x_0
x0为中心
λ
\lambda
λ为带宽内的样本点个数。不难发现这个方法得到的密度函数是一个阶梯函数(直方图),一种更加smooth的方法是基于核函数:
f
^
X
(
x
0
)
=
1
N
∑
i
=
1
N
K
λ
(
x
0
,
x
i
)
,
\hat{f}_X(x_0)=\frac{1}{N}\sum_{i=1}^NK_\lambda(x_0,x_i),
f^X(x0)=N1i=1∑NKλ(x0,xi),
在这个场景下最常用的核函数Gaussian kernel,也就是:
f
^
X
(
x
)
=
1
N
∑
i
=
1
N
ϕ
λ
(
x
−
x
i
)
=
(
F
^
⋆
ϕ
λ
)
(
x
)
,
\begin{array}{rcl}\hat{f}_X(x)&=&\frac{1}{N}\sum_{i=1}^N\phi_\lambda(x-x_i)\\&=&(\hat{F}\star\phi_\lambda)(x),\end{array}
f^X(x)==N1∑i=1Nϕλ(x−xi)(F^⋆ϕλ)(x),
后一个式子表示经验分布函数和
ϕ
\phi
ϕ的卷积。并且这个形式很容易可以推广到高维:
f
^
X
(
x
0
)
=
1
N
(
2
λ
2
π
)
p
2
∑
i
=
1
N
e
−
1
2
(
∣
∣
x
i
−
x
0
∣
∣
/
λ
)
2
.
\hat{f}_X(x_0)=\frac1{N(2\lambda^2\pi)^{\frac p2}}\sum_{i=1}^Ne^{-\frac12(||x_i-x_0||/\lambda)^2}.
f^X(x0)=N(2λ2π)2p1i=1∑Ne−21(∣∣xi−x0∣∣/λ)2.
Kernel Density Classification
使用核密度估计的方法得到核密度估计函数,从而基于贝叶斯的思路可以很容易得到given标签y的条件密度X|Y,再结合先验概率得到Y|X的后验概率:
Pr
^
(
G
=
j
∣
X
=
x
0
)
=
π
^
j
f
^
j
(
x
0
)
∑
k
=
1
J
π
^
k
f
^
k
(
x
0
)
.
\hat{\Pr}(G=j|X=x_0)=\frac{\hat{\pi}_j\hat{f}_j(x_0)}{\sum_{k=1}^J\hat{\pi}_k\hat{f}_k(x_0)}.
Pr^(G=j∣X=x0)=∑k=1Jπ^kf^k(x0)π^jf^j(x0).
但以上分类器需要建立在核密度估计效果较好的前提下,如果在边界点或者sparse的地方,那么核密度估计的方差就会很大,影响分类的效果。事实上,如果分类才是最终的目的,也不太建议使用这个方法来进行。
The Naive Bayes Classifier
朴素贝叶斯方法在处理高维、多种数据类型的分类问题上具有较大的优势,它做出假设:X的每一维都是独立的。这样的假设使得在高维下可以分别估计每一维的核密度再相乘即得到最终的核密度,以简化计算,并且可以包容X是定类的数据类型:
f
j
(
X
)
=
∏
k
=
1
p
f
j
k
(
X
k
)
.
f_j(X)=\prod_{k=1}^pf_{jk}(X_k).
fj(X)=k=1∏pfjk(Xk).
虽然每一维独立的假设看起来不太合适,但在实践中,它的表现比很多复杂的分类器效果都要好很多。同样的,可以使用对数似然比最大来判别所属类别:
log
Pr
(
G
=
ℓ
∣
X
)
Pr
(
G
=
J
∣
X
)
=
log
π
ℓ
f
ℓ
(
X
)
π
J
f
J
(
X
)
=
log
π
ℓ
∏
k
=
1
p
f
ℓ
k
(
X
k
)
π
J
∏
k
=
1
p
f
J
k
(
X
k
)
=
log
π
ℓ
π
J
+
∑
k
=
1
p
log
f
ℓ
k
(
X
k
)
f
J
k
(
X
k
)
=
α
ℓ
+
∑
k
=
1
p
g
ℓ
k
(
X
k
)
.
\begin{aligned} \log\frac{\operatorname*{Pr}(G=\ell|X)}{\operatorname*{Pr}(G=J|X)}& =\log\frac{\pi_\ell f_\ell(X)}{\pi_Jf_J(X)} \\ &=\log\frac{\pi_\ell\prod_{k=1}^pf_{\ell k}(X_k)}{\pi_J\prod_{k=1}^pf_{Jk}(X_k)} \\ &=\log\frac{\pi_\ell}{\pi_J}+\sum_{k=1}^p\log\frac{f_{\ell k}(X_k)}{f_{Jk}(X_k)} \\ &=\alpha_\ell+\sum_{k=1}^pg_{\ell k}(X_k). \end{aligned}
logPr(G=J∣X)Pr(G=ℓ∣X)=logπJfJ(X)πℓfℓ(X)=logπJ∏k=1pfJk(Xk)πℓ∏k=1pfℓk(Xk)=logπJπℓ+k=1∑plogfJk(Xk)fℓk(Xk)=αℓ+k=1∑pgℓk(Xk).
可以看到,它具有一般线性模型的形式,但是和局部逻辑回归不一样的是,局部逻辑回归虽然也具有这样的似然比函数形式,但是对于X|Y的分布是没有任何限制的,但这里使用核密度估计的方法进行了限制。朴素贝叶斯实际上是针对
P
(
X
∣
Y
=
i
)
P(X|Y=i)
P(X∣Y=i)进行建模的,由此得到
P
(
Y
∣
X
)
P(Y|X)
P(Y∣X);而逻辑回归则是直接针对
P
(
Y
∣
X
)
P(Y|X)
P(Y∣X)建模。
此外,朴素贝叶斯和general additive model(GAM)的区别在于,朴素贝叶斯假设每一维X都是独立的,但是GAM只是假设了相加的形式,这会导致在参数估计的时候出现查别。GAM往往需要使用backfitting或者是牛顿迭代来拟合,但是朴素贝叶斯基于独立性假设,使得可以对每一个维度单独使用核密度估计等方法来得到,不需要使用复杂的迭代算法。
Radial Basis Functions and Kernels
回顾第五章使用到的basis function的方法,由于knots的引入,实际上也具有一定的局部性,接下来希望可以将核函数的方法与basis function进行结合,使用核函数作为basis function,也就是radial basis functions:
f
(
x
)
=
∑
j
=
1
M
K
λ
j
(
ξ
j
,
x
)
β
j
=
∑
j
=
1
M
D
(
∣
∣
x
−
ξ
j
∣
∣
λ
j
)
β
j
,
\begin{array}{rcl}f(x)&=&\sum_{j=1}^MK_{\lambda_j}\left(\xi_j,x\right)\beta_j\\&=&\sum_{j=1}^MD\left(\frac{||x-\xi_j||}{\lambda_j}\right)\beta_j,\end{array}
f(x)==∑j=1MKλj(ξj,x)βj∑j=1MD(λj∣∣x−ξj∣∣)βj,
在这种应用场景下常用的核函数仍然是高斯核。
接下来介绍参数
(
λ
i
,
ξ
i
,
β
i
)
(\lambda_i,\xi_i,\beta_i)
(λi,ξi,βi)的估计方法,常用的一种还是最小二乘:
min
{
λ
j
,
ξ
j
,
β
j
}
1
M
∑
i
=
1
N
(
y
i
−
β
0
−
∑
j
=
1
M
β
j
exp
{
−
(
x
i
−
ξ
j
)
T
(
x
i
−
ξ
j
)
λ
j
2
}
)
2
.
\min_{\{\lambda_j,\xi_j,\beta_j\}_1^M}\sum_{i=1}^N\left(y_i-\beta_0-\sum_{j=1}^M\beta_j\exp\left\{-\frac{(x_i-\xi_j)^T(x_i-\xi_j)}{\lambda_j^2}\right\}\right)^2.
{λj,ξj,βj}1Mmini=1∑N(yi−β0−j=1∑Mβjexp{−λj2(xi−ξj)T(xi−ξj)})2.
-
如果 ξ i , λ i \xi_i,\lambda_i ξi,λi已经given,那么估计 β 0 , β i \beta_0,\beta_i β0,βi就是OLSE
-
估计 ξ i , λ i \xi_i,\lambda_i ξi,λi:只看X的分布特征,不考虑Y
例如可以使用 ξ i , λ i \xi_i,\lambda_i ξi,λi作为混合高斯的均值和标准差,从而使用EM算法进行估计;又例如可以使用聚类算法,从而得到每一类的均值和标准差估计值作为参数估计值;甚至可以直接把 λ i = λ \lambda_i=\lambda λi=λ作为超参数,对于 ξ i \xi_i ξi作为使用聚类算法得到聚类的中心点
虽然将
λ
i
=
λ
\lambda_i=\lambda
λi=λ可以减少参数估计的数量,但是这样每一个核就没有在宽度上的异质性了,会产生低谷\ref{Gaussian V.S. Renormalized Gaussian};更常用的方法是使用Renormalize radial basis functions:
h
j
(
x
)
=
D
(
∣
∣
x
−
ξ
j
∣
∣
/
λ
)
∑
k
=
1
M
D
(
∣
∣
x
−
ξ
k
∣
∣
/
λ
)
,
h_j(x)=\frac{D(||x-\xi_j||/\lambda)}{\sum_{k=1}^MD(||x-\xi_k||/\lambda)},
hj(x)=∑k=1MD(∣∣x−ξk∣∣/λ)D(∣∣x−ξj∣∣/λ),
如果使用Renormalized radial functions作为basis functions,则不难发现Nadaraya–Watson kernel regression estimator是它的一种特殊情况:
f
^
(
x
0
)
=
∑
i
=
1
N
y
i
K
λ
(
x
0
,
x
i
)
∑
i
=
1
N
K
λ
(
x
0
,
x
i
)
=
∑
i
=
1
N
y
i
h
i
(
x
0
)
\begin{aligned} \hat{f}(x_{0}) =\sum_{i=1}^{N}y_{i}\frac{K_{\lambda}\left(x_{0},x_{i}\right)}{\sum_{i=1}^{N}K_{\lambda}(x_{0},x_{i})} \\ =\sum_{i=1}^{N}y_{i}h_{i}(x_{0}) \end{aligned}
f^(x0)=i=1∑Nyi∑i=1NKλ(x0,xi)Kλ(x0,xi)=i=1∑Nyihi(x0)
可以认为相当于
λ
^
i
=
λ
,
ξ
^
i
=
x
i
,
β
^
i
=
y
i
\hat \lambda_i=\lambda,\hat \xi_i=x_i,\hat \beta_i=y_i
λ^i=λ,ξ^i=xi,β^i=yi。
Radial basis functions搭建了核方法和局部回归之间的桥梁。
Mixture Models for Density Estimation and Classification
高斯混合模型密度函数形式:
f
(
x
)
=
∑
m
=
1
M
α
m
ϕ
(
x
;
μ
m
,
Σ
m
)
f(x)=\sum_{m=1}^M\alpha_m\phi(x;\mu_m,\Sigma_m)
f(x)=m=1∑Mαmϕ(x;μm,Σm)
其中
∑
m
=
1
M
α
m
=
1
\sum_{m=1}^M\alpha_m=1
∑m=1Mαm=1。这可以看成是核方法的一个推广:
-
如果 Σ m = σ m I \Sigma_m=\sigma_mI Σm=σmI,就相当于一个 λ m = σ m , ∑ m = 1 M α m = 1 \lambda_m=\sigma_m,\sum_{m=1}^M\alpha_m=1 λm=σm,∑m=1Mαm=1的radial basis expansion
-
如果 σ m = σ \sigma_m=\sigma σm=σ是固定的,M和N相等,则以上估计的结果会接近于 α ^ m = 1 / N , μ ^ m = x m \hat \alpha_m=1/N,\hat \mu_m=x_m α^m=1/N,μ^m=xm
使用混合模型对具有分类标签的总体进行核密度估计,虽然没有具体使用标签,但是仍然可以建模出两种总体的区别。