关于样条 - 经管之家

关于样条

zjj8211

3210

收藏 2008-10-29

请教各位高手,怎么样用R来实现用样条(比如B样条)来渐进回归函数

[em01]

扫码加我拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

全部回复

DM小菜鸟

2014-12-24 18:03:00

样条估计
如果函数在不同地方有不同的非线性度，或者有多个极值点，那么用多项式特别是低阶多项式来完成拟合是非常不合适的。一种解决办法是我们之前提到的近邻多项式（或者称局部多项式），另一种就是样条——用分段的低阶多项式逼近函数。
关于样条，常用的有两类，一类是多项式样条，另一类是光滑样条。

多项式样条多项式样条的样条基有很多，最为著名的是我们之前在函数逼近中提到的truncated power basis与B-spline basis。我们这里十分简要的介绍一下B样条，B样条基下的函数逼近可以写为：[ f(x)=beta_0+beta_1 x+cdots+beta_p x^p+sum_{j=1}^n beta_j B_j^p(x) ]其中[ B_i^p(x)=frac{x-c_i}{c_{i+p}-c_i}B_{i}^{p-1}(x)+frac{c_{i+p+1}-x}{c_{i+p+1}-c_{i+1}}B_{i+1}^{p-1}(x) ]上式中( B_i^0(x) =1 )当且仅当( c_i le x<c_{i+1} )否则取0.在R中splines包的函数bs()提供了B样条估计，其调用格式为：

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

对于参数df值得说明的是df=degree+（Knots个数），attr(,“knots”)会显示划分点，我们常用的3次B样条公式: df=k+3 (不含常数项)

essay data为例说明B样条的估计情况：

easy <- read.table("D:/R/data/easysmooth.dat", header = T)x <- easy$Xy <- easy$Ym.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")  #可以将看到，节点在不指定的情况下默认的是均匀样条，当然，我们可以根据散点图给#出节点的具体选择。  ## 25% 50% 75% ## -1.875 -0.250  1.375m.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.1AIC(m.bsp1)## [1] 727.4summary(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样条基本很接近真实函数了，summary(m.bsp)报告了各个系数的估计，带入( f(x) )的B样条基展开中即可得到一个显式的表达式。

光滑样条
虽然B样条已经很好了，但是理论与实践都表明直接用最小二乘去求解系数效果不好，容易过拟合。一个可能的改进是光滑样条。所谓的光滑样条，就是在求解最小二乘时给估计函数( f(x) )加上了一定的惩罚，这个有点类似压缩估计。我们这里采用最常用的光滑性惩罚，得到函数( f(x) )的估计( m(x) )满足如下的惩罚最小二乘：[ min sum_{i=1}^n (y_i-m(x_i))^2+lambda int [m''(x)]^2 dx ]在R的splines包中提供了函数smooth.spline来求解光滑样条

easy <- read.table("D:/R/data/easysmooth.dat", header = T)x <- easy$Xy <- easy$Ys.hat <- smooth.spline(x, y)## OUTPUTs.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 PLOTSs <- 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)

最后我们来讲一下怎么计算出( m(x) ),这里我们使用Reinsch algorithm。Step 1: 计算向量( Q'y ) .Step 2: 找到一个非0对角阵( R+lambda Q'Q ) 使得它可以进行Cholesky分解，有因子L，DStep 3: 解方程：( (R+lambda Q'Q)gamma=Q'y )Step 4: 得到估值( m=y-alpha Q gamma ).上面的Q与R可以表示为：

上面的t表示节点。我们不妨来算算essay data的例子：
easy <- read.table("D:/R/data/easysmooth.dat", header = T)x <- easy$Xy <- easy$Yn <- 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 Q[i + 1, i] = -1/h - 1/h[i + 1] Q[i + 2, i] = 1/h[i + 1]}for (i in 2:(n - 2)) { R[i, i] = 1/6 * (h + h[i + 1]) R[i - 1, i] = h/6 R[i, i - 1] = h/6}R[1, 1] = 1/6 * (h[1] + h[2])lambda <- 0.2A <- R + lambda * t(Q) %*% Qgamma <- solve(A, t(Q) %*% as.matrix(y))g <- as.matrix(y) - lambda * Q %*% gammas <- 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)
在惩罚系数为0.2的情况下，拟合还是不坏的，不是吗？至于为什么可以这样算，我们只要注意到( int [m^{''}(x)]dx=m^'(x_i)QR^{-1}Q^'m(x_i) ),估计的问题就与我们十分熟悉的lasso，岭回归十分相像了。

扫码加我拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

扫码加我 拉你入群

扫码加我 拉你入群

分享

扫码加好友，拉您进群

扫码加我拉你入群

扫码加我拉你入群