1.前言
前面介绍spline的基本性质和最小二乘回归，这里进入正题，主要总结Spline Regression的内容，主要的复现参考论文中的P-spline,针对的是Cardinal B-spline，就是分割节点是等间隔的
2.基本理论
高阶基函数
B
i
,
k
(
x
)
{{B}_{i,k}}(x)
Bi,k(x)的递推公式
B
i
,
k
(
x
)
:
=
x
−
t
i
t
i
+
k
−
t
i
B
i
,
k
−
1
(
x
)
+
t
i
+
k
+
1
−
x
t
i
+
k
+
1
−
t
i
+
1
B
i
+
1
,
k
−
1
(
x
)
.
{{B}_{i,k}}(x):=\frac{x-{{t}_{i}}}{{{t}_{i+k}}-{{t}_{i}}}{{B}_{i,k-1}}(x)+\frac{{{t}_{i+k+1}}-x}{{{t}_{i+k+1}}-{{t}_{i+1}}}{{B}_{i+1,k-1}}(x).
Bi,k(x):=ti+k−tix−tiBi,k−1(x)+ti+k+1−ti+1ti+k+1−xBi+1,k−1(x).
B
i
,
k
(
x
)
{{B}_{i,k}}(x)
Bi,k(x)的一阶导数
d
B
i
,
k
(
x
)
d
x
=
k
(
B
i
,
k
−
1
(
x
)
t
i
+
k
−
t
i
−
B
i
+
1
,
k
−
1
(
x
)
t
i
+
k
+
1
−
t
i
+
1
)
\frac{d{{B}_{i,k}}(x)}{dx}=k\left( \frac{{{B}_{i,k-1}}(x)}{{{t}_{i+k}}-{{t}_{i}}}-\frac{{{B}_{i+1,k-1}}(x)}{{{t}_{i+k+1}}-{{t}_{i+1}}} \right)
dxdBi,k(x)=k(ti+k−tiBi,k−1(x)−ti+k+1−ti+1Bi+1,k−1(x))
在Cardinal B-spline时，假设间隔为
h
h
h,则有
d
B
i
,
k
(
x
)
d
x
=
1
h
(
B
i
,
k
−
1
(
x
)
−
B
i
+
1
,
k
−
1
(
x
)
)
\frac{d{{B}_{i,k}}(x)}{dx}=\frac{1}{h}\left( {B}_{i,k-1}(x)-{{B}_{i+1,k-1}}(x)\right)
dxdBi,k(x)=h1(Bi,k−1(x)−Bi+1,k−1(x))
B-spline的加权最小二乘法
α
^
=
a
r
g
m
i
n
α
∑
all
x
{
W
(
x
)
[
y
(
x
)
−
∑
i
α
i
B
i
,
k
,
t
(
x
)
]
}
2
\hat{ \mathbf{\alpha}}=\underset{\alpha}{argmin} \underset{\text{all}\,x}{ \sum }{{\left\{ W(x)\left[ y(x)-\underset{i}{ \sum }\,{{\alpha }_{i}}{{B}_{i,k,t}}(x) \right] \right\}}^{2}}
α^=αargminallx∑{W(x)[y(x)−i∑αiBi,k,t(x)]}2
W
(
x
)
W(x)
W(x)是一个权值系数，这里取单位矩阵
I
I
I,
α
i
{\alpha }_{i}
αi为待定的最小二乘系数
3惩罚项
3.1惩罚项介绍
为了改进这些最小二乘系数呢，在最小二乘的目标式子中，加入惩罚项。不同的惩罚项有不同的效果。
从前面的最小二乘法我们可以得到B-spline的曲线拟合
y
^
\hat{y}
y^
y
^
=
∑
α
i
B
i
,
k
(
x
)
\hat{y} = { \sum }\,{{\alpha }_{i}}{{B}_{i,k}}(x)
y^=∑αiBi,k(x)
惩罚项可以从
y
^
\hat{y}
y^的导数连续程度，以及
α
\alpha
α的连续程度两个来入手
y
^
\hat{y}
y^的一阶导数
d
y
^
d
x
=
d
d
x
∑
i
α
i
B
i
,
k
=
∑
i
α
i
B
i
,
k
′
=
−
1
h
∑
j
Δ
α
j
+
1
B
j
(
x
;
q
−
1
)
\frac{d\hat{y}}{dx}=\frac{d}{dx}\underset{i}{ \sum }\,{{\alpha }_{i}}{{B}_{i,k}}=\underset{i}{ \sum }\,{{\alpha }_{i}}{{B}_{i,k}^{\prime}}=-\frac{1}{h}\sum_j{\Delta \alpha_{j+1}}B_{j}(x;q-1)
dxdy^=dxdi∑αiBi,k=i∑αiBi,k′=−h1j∑Δαj+1Bj(x;q−1)
y
^
\hat{y}
y^的二阶导数
∑ j a j B j ′ ′ ( x ; q ) = h − 2 ∑ j Δ 2 a j B j ( x ; q − 2 ) Δ 2 a j = Δ Δ a j = a j − 2 a j − 1 + a j − 2 \sum_{j} a_{j} B_{j}^{\prime \prime}(x ; q)=h^{-2}\sum_{j} \Delta^{2} a_{j} B_{j}(x ; q-2)\\ \Delta^{2} a_{j}=\Delta \Delta a_{j}=a_{j}-2 a_{j-1}+a_{j-2} j∑ajBj′′(x;q)=h−2j∑Δ2ajBj(x;q−2)Δ2aj=ΔΔaj=aj−2aj−1+aj−2
3.2惩罚项
3.2.1 二阶导数
S = ∑ i = 1 m { y i − ∑ j = 1 n a j B j ( x i ) } 2 + λ ∫ x min x max { ∑ j = 1 n a j B j ′ ′ ( x ) } 2 d x \begin{aligned} S=& \sum_{i=1}^{m}\left\{y_{i}-\sum_{j=1}^{n} a_{j} B_{j}\left(x_{i}\right)\right\}^{2} +\lambda \int_{x_{\min }}^{x_{\max }}\left\{\sum_{j=1}^{n} a_{j} B_{j}^{\prime \prime}(x)\right\}^{2} d x \end{aligned} S=i=1∑m{yi−j=1∑najBj(xi)}2+λ∫xminxmax{j=1∑najBj′′(x)}2dx
对二阶导的积分，其实只是对二导求一个和。最终会以离散的形式计算，而不是以积分的形式。这个和越大，代表曲线越不光滑。
举例来说，以3阶B-spline的二阶导数作为惩罚项
h
2
P
=
λ
∫
x
min
x
max
{
∑
j
a
j
B
j
′
′
(
x
;
3
)
}
2
d
x
=
λ
∫
x
min
x
max
{
∑
j
Δ
2
a
j
B
j
(
x
;
1
)
}
2
d
x
=
λ
∫
x
min
x
max
∑
j
∑
k
Δ
2
a
j
Δ
2
a
k
⋅
B
j
(
x
;
1
)
B
k
(
x
;
1
)
d
x
h^{2} P=\lambda \int_{x_{\min }}^{x_{\max }}\left\{\sum_{j} a_{j} B_{j}^{\prime \prime}(x ; 3)\right\}^{2} d x= \lambda \int_{x_{\min }}^{x_{\max }}\left\{\sum_{j} \Delta^{2} a_{j} B_{j}(x ; 1)\right\}^{2} d x\\= \lambda \int_{x_{\min }}^{x_{\max }} \sum_{j} \sum_{k} \Delta^{2} a_{j} \Delta^{2} a_{k} \cdot B_{j}(x ; 1) B_{k}(x ; 1) d x
h2P=λ∫xminxmax{j∑ajBj′′(x;3)}2dx=λ∫xminxmax{j∑Δ2ajBj(x;1)}2dx=λ∫xminxmaxj∑k∑Δ2ajΔ2ak⋅Bj(x;1)Bk(x;1)dx
对于一阶B-spline的基函数
B
j
(
x
;
1
)
B_{j}(x ; 1)
Bj(x;1),它的作用范围为只在
[
t
j
,
t
j
+
1
)
[t_j,t_{j+1})
[tj,tj+1),那么，
k
k
k只有在
k
∈
[
j
−
1
,
j
+
1
]
k \in [j-1,j+1]
k∈[j−1,j+1]，两者才存在overlap，其余皆为0.上式可以继续化简，但是最终的结果依旧比较复杂。这就是文献1，为什么要提出了系数高阶导的方法
3.2.2 系数高阶导
S
=
∑
i
=
1
m
{
y
i
−
∑
j
=
1
n
a
j
B
j
(
x
i
)
}
2
+
λ
∑
j
=
k
+
1
n
(
Δ
k
a
j
)
2
S=\sum_{i=1}^{m}\left\{y_{i}-\sum_{j=1}^{n} a_{j} B_{j}\left(x_{i}\right)\right\}^{2}+\lambda \sum_{j=k+1}^{n}\left(\Delta^{k} a_{j}\right)^{2}
S=i=1∑m{yi−j=1∑najBj(xi)}2+λj=k+1∑n(Δkaj)2
写成矩阵的矩阵形式如下
B
T
y
=
(
B
T
B
+
λ
D
k
T
D
k
)
a
B^{T} y=\left(B^{T} B+\lambda D_{k}^{T} D_{k}\right) a
BTy=(BTB+λDkTDk)a
B基函数矩阵，每一列代表一个基函数产生的序列，
D
k
D_k
Dk代表了系数
a
\mathbf a
a的k阶差分矩阵，且有
Δ
a
k
=
D
k
a
\Delta a^k=D_{k}a
Δak=Dka.在计算
Δ
a
k
\Delta a^k
Δak 时，这个看似不好表达的惩罚项目，亏了优美的线性代数，这项还是容易表达的。只要对单位矩阵
I
k
I_k
Ik做k次按行差分即可，matlab和R语言都有内置diff函数可以直接用。需要要注意的是，差分项随着差分阶数的增加而减少。如果没有
D
k
D_k
Dk，这个方程就是最普通的最小二乘回归。
在论文采用的是系数的2阶导数，到了这里，我们已经可以可以利用
λ
\lambda
λ控制拟合的平滑程度。如何选择一个合适的
λ
\lambda
λ，需要有了评价指标。这篇文章选择了
A
k
a
i
k
e
i
n
f
o
r
m
a
t
i
o
n
c
r
i
t
e
r
i
o
n
(
A
I
C
)
Akaike\ information \ criterion(AIC)
Akaike information criterion(AIC)指标，属于信息论里的概念，有许多种定义，目测还不等价，我引用了wiki中的定义。
A
I
C
=
2
k
−
2
L
L
=
−
l
n
(
(RSS/n)
)
⇒
A
I
C
=
2
k
+
l
n
(
R
S
S
/
n
)
A I C=2 k-2 L\\ L= -ln (\sqrt{\operatorname{(RSS/n)}} ) \\ \Rightarrow AIC = 2k+ln(RSS/n)\\
AIC=2k−2LL=−ln((RSS/n))⇒AIC=2k+ln(RSS/n)
k
k
k是参数数量，书中指的是有效参数数量，计算比想象得复杂得多，
R
S
S
RSS
RSS残差，
n
n
n代表样本数量。意义很好理解，显然AIC越小越好。
论文中的定义更加复杂，不知道是不是扯上信息论，都比较复杂
书中最终采用是GCV而不是直接使用AIC，其计算过程
H
=
B
(
B
T
W
~
B
+
λ
D
k
T
D
k
)
−
1
B
T
W
~
H=B\left(B^{T} \tilde{W} B+\lambda D_{k}^{T} D_{k}\right)^{-1} B^{T} \tilde{W}
H=B(BTW~B+λDkTDk)−1BTW~
W
~
\tilde{W}
W~为估计的加权系数
GCV计算过程
GCV
(
λ
)
=
∑
i
=
1
m
(
y
i
−
y
^
i
)
2
(
m
−
∑
i
=
1
m
h
i
i
)
2
\operatorname{GCV}(\lambda)=\sum_{i=1}^{m} \frac{\left(y_{i}-\hat{y}_{i}\right)^{2}}{\left(m-\sum_{i=1}^{m} h_{i i}\right)^{2}}
GCV(λ)=i=1∑m(m−∑i=1mhii)2(yi−y^i)2
这公式分子部分代表了残差，分母代表了误差的方差估计。论文假定误差是正态分布的。这公式中，
h
i
i
h_{i i}
hii体现出有效参数数量k.因为
λ
\lambda
λ比较大的时候，系数会被压缩，它的有效长度会减少。为什么用这个来表征k呢，这个细节我没弄明白。
T
(
λ
)
=
∑
i
=
1
m
h
i
i
=
tr
(
H
)
=
tr
{
(
I
+
λ
L
)
−
1
}
=
∑
j
=
1
n
1
1
+
λ
γ
j
T(\lambda)=\sum_{i=1}^{m} h_{i i}=\operatorname{tr}(H)=\operatorname{tr}\left\{(I+\lambda L)^{-1}\right\}=\sum_{j=1}^{n} \frac{1}{1+\lambda \gamma_{j}}
T(λ)=i=1∑mhii=tr(H)=tr{(I+λL)−1}=j=1∑n1+λγj1
γ
j
\gamma_{j}
γj是
L
L
L矩阵的特征向量，
L
L
L矩阵是
H
H
H的一部分，具体就不展开了
我感觉作者总想使算法变得更加generalized，所以无形中，感觉简单的部分变得复杂起来。
4.复现代码
论文中，附加了S-PLUS和matlab代码。所以，也谈不上说复现，他用到的测试数据我也没有去找，自己生成一个数据就完事了。
我就用他的S-PLUS代码，稍作修改一下，在R环境里面就可以跑了
作者的代码比 文章要简单太多，也暴露自己基础不扎实的问题
bspline<-function(x,xl,xr,ndx,bdeg){
dx<-(xr-xl)/ndx
knots<-seq(xl-bdeg*dx,xr+bdeg*dx,by=dx)
B<-spline.des(knots,x,bdeg+1,0*x)$design
B
}
pord <-2
n <- 100
x <- seq(0, 1, length=n)
knots <-seq(0, 1, length=10)
bdeg <-3
lambda <-100
dgp <- sin(2*pi*x)
y <- dgp + rnorm(n, sd=.2)
B <- spline.des(knots,x, bdeg,0,outer.ok = TRUE)
B <- B$design
D <- diag(ncol(B))
yhat<-NULL
diff(D,differences = pord)
gcv <- 1:10
lambdaArray <- 1:10
for(i in 1:10){
lambda = (2^i-2)/50
lambdaArray[i] = lambda
a <- solve(t(B) %*% B + lambda * t(D) %*% D,t(B) %*% y)
yhat<-cbind(yhat,B %*% a)
s <- sum((y - yhat)^2)
Q <- solve(t(B) %*% B + lambda * t(D) %*% D)
# matrix inversion
t <- sum(diag(Q %*% (t(B) %*% B)))
gcv[i] <- s / (nrow(B) - t)^2
}
par(mfrow = c(3,1))
matplot(x, B, type="l", lwd=2)
plot(x, y, cex=.25, col="grey")
lines(x, dgp, lty=2)
lines(x, yhat[,1], col="red",lwd=0.5)
lines(x, yhat[,5], col="blue",lwd=0.5)
lines(x, yhat[,9], col="green", lwd=0.5)
legend("topright",legend=c("s+n","s","lb=0","lb=0.6","lb=10.2"), col=c("grey","black","red","blue","green"), lty=1,lwd=2)
plot(lambdaArray, gcv,type = "b")
第一幅图是基函数，第二幅图是不同
λ
\lambda
λ的拟合情况情况，
λ
\lambda
λ越大，曲线越光整，简直到了平整，绿色的线到了往水平直线方向趋近了。第三幅图可以看到gcv曲线是单调上升的。也就是说这个数据其实并不需要做惩罚，本身已经达到很好的拟合效果。
5.小结
这篇论文写得不错，引用量很高。spline的回归和惩罚讲得比较清楚，特别是背后的原理部分，有许多值得借鉴。作者在写作的时候，感觉还是个小年轻，措辞能力上感觉还有待改进。但是可以看到作者的理论功力还是不错，有许多细节值得深挖下去，论文后面带的评价也比较高。由于研究的方向和作者的不太相同，不想花太多时间。不过花了这么几天功夫去学习，还是有一些收获。
References
- Eilers, P.H.C. and Marx, B.D. (1996). Flexible smoothing with B-splines and penalties (with comments and rejoinder). Statistical Science 11(2): 89-121.
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
