ESL-chapter5 Smoothing Splines

先来看公式5.9

lamda称之为光滑系数。第一项衡量拟合函数的近似程度，第二项惩罚不光滑的拟合函数。一般来说，一阶导数衡量斜率的变化，二阶导数衡量斜率的程度~曲线变动是否剧烈。对二阶导数的积分视作为斜率变化程度的累加和。

如果lamda=infinity，则函数f()的二阶微分为0，说明f（）为线性函数。

如果lamda=0，则f()完全拟合各个样本点，形成剧烈变动的拟合函数（过拟合）。

满足5.9的函数回归基函数需要满足如下几个条件：

1，knots需要落在样本点上。（可以计算一阶二阶导数）

2，一阶二阶导数存在。

3，样本点以外的函数要是线性的（使得第二项为0）。

因此自然三次样条是一个很显然的解，其knots定义在每个样本点上。但是由于lamda的存在，该自然三次样条是一个朝着线性函数收缩的版本。

下面来看图5.6如何绘制。

代码如下

library(ElemStatLearn) 

males = bone$gender == "male" 

females = bone$gender == "female"
boneMaleSmooth = smooth.spline( bone[males,"age"], bone[males,"spnbmd"], df=12 )#前一个参数是自变量，后一个是因变量。df相当于选择了lamda，定义在下面介绍。
boneFemaleSmooth = smooth.spline( bone[females,"age"], bone[females,"spnbmd"], df=12 )
plot(boneMaleSmooth, ylim=c(-0.05,0.20), col="blue", type="l", xlab="Age", ylab="spnbmd")
points(bone[males,c(2,4)], col="blue", pch=20)
lines(boneFemaleSmooth, ylim=c(-0.05,0.20), col="red")
points(bone[females,c(2,4)], col="red", pch=20)

下面看5.4.1节

比较重要的概念是有效自由度，定义为光滑矩阵的迹的和。直接决定了拟合曲线的光滑程度。来看图5.7的代码

LAozone = read.table("LAozone.data",sep=",",head=T)
attach(LAozone)
train=sample(330,128)#采样
LAozoneSmooth.samller = smooth.spline( LAozone[train,"dpg"], LAozone[train,"ozone"], df=11)# 有效自由度为11
LAozoneSmooth.bigger = smooth.spline( LAozone[train,"dpg"], LAozone[train,"ozone"], df=5)# 有效自由度为5
plot(LAozoneSmooth.samller, ylim=range(LAozone$ozone), col="green", type="l", xlab="Daggot Pressure Gradient", ylab="Ozone Concentration")
points( LAozone[train,c("dpg","ozone")], col="black", pch=20)
lines(LAozoneSmooth.bigger, ylim=range(LAozone$ozone), col="red", type="l")

从图中可以看出较大自由度的曲线波动较大。

把光滑矩阵Slamda进行特征值分解的目的是看Slamda是如何作用于y,从而影响拟合函数的。光滑矩阵和y没有关系，只和x（一维的）和自由度有关。

下面，贴出如何求光滑矩阵的代码

smooth.matrix = function(x,df){
## return the smoother matrix with knots x and degree of freedom = df
##该函数要求x的值是唯一的，并且是排序之后的。
 n = length(x);
 A = matrix(0, n, n);
 for(i in 1:n){
       y = rep(0, n); y[i]=1;
       yi = smooth.spline(x, y, df=df)$y;
       A[,i]= yi;
}
 (A+t(A))/2;
}

S.LAozone.5 <- smooth.matrix(sort(unique(LAozone[train,"dpg"])),df=5)#求自由度为5的大气污染数据的光滑矩阵
S.LAozone.11 <- smooth.matrix(sort(unique(LAozone[train,"dpg"])),df=11)
df=0; for(i in 1:81) df=df+S.LAozone.5[i,i];
df#验证一下该计算是对的。df=5.000571
然后根据特征值分解

 eigen.S.LAozone.5 = eigen(S.LAozone.5);
eigen.S.LAozone.11 = eigen(S.LAozone.11);
## plot the eigen values
## 画出书上图5.7的左下角的图
plot(eigen.S.LAozone.5$va[1:25], pch=1, col="red",ylab="Eigenvalues");
points(eigen.S.LAozone.11$va[1:25], pch=4, col="green");
lines(c(0,25), c(1, 1), col="black");

特征向量的图，我画的和书上不太一样，不知道为什么，代码就不贴了。

5.5 选择光滑参数

先上一个画图5.8的代码

myimage=function(x, mycol=heat.colors(12)){
  ## plot the matrix x as an image
  nr=dim(x)[1];
  tmp=x;
  for(i in 1:nr) tmp[i,]=x[(nr-i+1),];
  image(t(tmp), col=mycol)
}
myimage(S.LAozone.5 , mycol=terrain.colors(12))

光滑参数直接决定了回归基函数的次数，knot的个数和位置。因此省去了knot的选择问题。但是需要选择lamda,由于lamda和曲线的光滑程度是单调正相关的，因此可以用交叉检验来比较选择lamda.再配合可加性模型，就可以进行解决多变量的拟合问题。

5.5.2 节是方差偏倚分解。

首先，我们根据sin函数生成拟合数据。代码如下：

x <- seq(from =0, to=1, by=0.01)
f <- sin(12*(x+0.2))/(x+0.2)
y <- f+rnorm(100,0,1)#按照书上的公式生成数据。

其中，黄色曲线是真实数据，蓝色是df=9,绿色是df=5,红色是df=15，从中我们可以看出df=9和真实数据最接近。

生成置信带的代码如下 df=15

s.fit.15 <- smooth.matrix(x,df=15)
cov.f.15 <- s.fit.15%*%t(s.fit.15)
se.15 <- sqrt(diag(cov.f.15))
pred.15=predict(y.fit.15,newdata=x,se=T)
plot(x,y)
lines(x,pred.15$y,col="green",lwd=4)
lines(x,pred.15$y-2*se.15,col="green",lwd=2,lty="dashed")
lines(x,pred.15$y+2*se.15,col="green",lwd=2,lty="dashed")

从书上的图可以看出，df=5时，bias较大，置信带较窄。而df=15时，bias较小，但是置信带较宽，变化很大。这就是偏倚方差分解。df=9时，较为折中。