R語(yǔ)言與函數(shù)估計(jì)學(xué)習(xí)筆記（樣條方法）

樣條估計(jì)

如果函數(shù)在不同地方有不同的非線性度，或者有多個(gè)極值點(diǎn)，那么用多項(xiàng)式特別是低階多項(xiàng)式來(lái)完成擬合是非常不合適的。一種解決辦法是我們之前提到的近鄰多項(xiàng)式（或者稱局部多項(xiàng)式），另一種就是樣條——用分段的低階多項(xiàng)式逼近函數(shù)。
關(guān)于樣條，常用的有兩類，一類是多項(xiàng)式樣條，另一類是光滑樣條。

多項(xiàng)式樣條

多項(xiàng)式樣條的樣條基有很多，最為著名的是我們之前在函數(shù)逼近中提到的truncated power basis與B-spline basis。我們這里十分簡(jiǎn)要的介紹一下B樣條，B樣條基下的函數(shù)逼近可以寫為：

其中

上式中否則取0.在R中splines包的函數(shù)bs()提供了B樣條估計(jì)，其調(diào)用格式為：

bs(x, df = NULL, knots = NULL, degree = 3, intercept = FALSE, Boundary.knots = range(x))

對(duì)于參數(shù)df值得說(shuō)明的是df=degree+（Knots個(gè)數(shù)），attr(,“knots”)會(huì)顯示劃分點(diǎn)，我們常用的3次B樣條公式: df=k+3 (不含常數(shù)項(xiàng))
我們以前面提到的essay data為例說(shuō)明B樣條的估計(jì)情況：

easy <- read.table("D:/R/data/easysmooth.dat", header = T)
x <- easy$X
y <- easy$Y
m.bsp <- lm(y ~ bs(x, df = 6))


s = function(x) {
    (x^3) * sin((x + 3.4)/2)
}
x.plot = seq(min(x), max(x), length.out = 1000)
y.plot = s(x.plot)
plot(x, y, xlab = "Predictor", ylab = "Response")
lines(x.plot, y.plot, lty = 1, col = 1)
lines(x, fitted(m.bsp), lty = 2, col = 2)

attr(bs(x, df = 6), "knots")  #可以將看到，節(jié)點(diǎn)在不指定的情況下默認(rèn)的是均勻樣條，當(dāng)然，我們可以根據(jù)散點(diǎn)圖給#出節(jié)點(diǎn)的具體選擇。

##    25%    50%    75%
## -1.875 -0.250  1.375


m.bsp1 <- lm(y ~ bs(x, df = 6, knots = c(-2.5, -1, 2)))
lines(x, fitted(m.bsp1), lty = 3, col = 3)

AIC(m.bsp)

## [1] 718.1

AIC(m.bsp1)

## [1] 727.4


summary(m.bsp)

##
## Call:
## lm(formula = y ~ bs(x, df = 6))
##
## Residuals:
##    Min     1Q Median     3Q    Max
## -3.790 -0.911 -0.065  0.892  4.445
##
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       1.816      0.622    2.92   0.0039 **
## bs(x, df = 6)1  -10.552      1.161   -9.09  < 2e-16 ***
## bs(x, df = 6)2   -7.127      0.755   -9.44  < 2e-16 ***
## bs(x, df = 6)3    0.813      0.926    0.88   0.3808    
## bs(x, df = 6)4   -4.056      0.859   -4.72  4.5e-06 ***
## bs(x, df = 6)5    5.781      0.967    5.98  1.1e-08 ***
## bs(x, df = 6)6   -3.505      0.865   -4.05  7.4e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.42 on 193 degrees of freedom
## Multiple R-squared:  0.824,  Adjusted R-squared:  0.819
## F-statistic:  151 on 6 and 193 DF,  p-value: <2e-16

可以看到B樣條基本很接近真實(shí)函數(shù)了，summary(m.bsp)報(bào)告了各個(gè)系數(shù)的估計(jì)，帶入f(x)的B樣條基展開(kāi)中即可得到一個(gè)顯式的表達(dá)式。

光滑樣條

雖然B樣條已經(jīng)很好了，但是理論與實(shí)踐都表明直接用最小二乘去求解系數(shù)效果不好，容易過(guò)擬合。一個(gè)可能的改進(jìn)是光滑樣條。所謂的光滑樣條，就是在求解最小二乘時(shí)給估計(jì)函數(shù)f(x)加上了一定的懲罰，這個(gè)有點(diǎn)類似壓縮估計(jì)。我們這里采用最常用的光滑性懲罰，得到函數(shù)f(x)的估計(jì)m(x)滿足如下的懲罰最小二乘：

在R的splines包中提供了函數(shù)smooth.spline來(lái)求解光滑樣條

easy <- read.table("D:/R/data/easysmooth.dat", header = T)
x <- easy$X
y <- easy$Y
s.hat <- smooth.spline(x, y)

## OUTPUT
s.hat

## Call:
## smooth.spline(x = x, y = y)
## 
## Smoothing Parameter  spar= 0.7251  lambda= 0.0002543 (12 iterations)
## Equivalent Degrees of Freedom (Df): 11.56
## Penalized Criterion: 380.9
## GCV: 2.145

## OUTPUT PLOTS
s <- function(x) {
    (x^3) * sin((x + 3.4)/2)
}
x.plot = seq(min(x), max(x), length.out = 1000)
y.plot = s(x.plot)
plot(x, y, xlab = "Predictor", ylab = "Response")
lines(x.plot, y.plot, lty = 1, col = 1)
lines(s.hat, lty = 2, col = 2)

最后我們來(lái)講一下怎么計(jì)算出m(x),這里我們使用Reinsch algorithm。Step 1: 計(jì)算向量Q′y.Step 2: 找到一個(gè)非0對(duì)角陣R+λQ′Q使得它可以進(jìn)行Cholesky分解，有因子L，DStep 3: 解方程：(R+λQ′Q)γ=Q′yStep 4: 得到估值m=y?αQγ.上面的Q與R可以表示為：

上面的t表示節(jié)點(diǎn)。我們不妨來(lái)算算essay data的例子：

easy <- read.table("D:/R/data/easysmooth.dat", header = T)
x <- easy$X
y <- easy$Y

n <- length(y)
knots <- seq(min(x), max(x), length = n + 1)
h <- knots[-1] - knots[-n]
Q <- matrix(0, n, n - 2)
R <- matrix(0, n - 2, n - 2)
for (i in 1:(n - 2)) {
    Q[i, i] = 1/h[i]
    Q[i + 1, i] = -1/h[i] - 1/h[i + 1]
    Q[i + 2, i] = 1/h[i + 1]
}
for (i in 2:(n - 2)) {
    R[i, i] = 1/6 * (h[i] + h[i + 1])
    R[i - 1, i] = h[i]/6
    R[i, i - 1] = h[i]/6
}
R[1, 1] = 1/6 * (h[1] + h[2])
lambda <- 0.2
A <- R + lambda * t(Q) %*% Q
gamma <- solve(A, t(Q) %*% as.matrix(y))

g <- as.matrix(y) - lambda * Q %*% gamma

s <- function(x) {
    (x^3) * sin((x + 3.4)/2)
}
x.plot <- seq(min(x), max(x), length.out = 1000)
y.plot <- s(x.plot)
plot(x, y, xlab = "Predictor", ylab = "Response")
lines(x.plot, y.plot, lty = 1, col = 1)
lines(x, g, lty = 2, col = 2)

在懲罰系數(shù)為0.2的情況下，擬合還是不壞的，不是嗎？至于為什么可以這樣算，我們只要注意到\int [m^{''}(x)]dx=m^'(x_i)QR^{-1}Q^'m(x_i),估計(jì)的問(wèn)題就與我們十分熟悉的lasso，嶺回歸十分相像了。