R语言与函数估计学习笔记(样条方法)
2014-05-17 11:32
543 查看
样条估计
如果函数在不同地方有不同的非线性度,或者有多个极值点,那么用多项式特别是低阶多项式来完成拟合是非常不合适的。一种解决办法是我们之前提到的近邻多项式(或者称局部多项式),另一种就是样条——用分段的低阶多项式逼近函数。关于样条,常用的有两类,一类是多项式样条,另一类是光滑样条。多项式样条
多项式样条的样条基有很多,最为著名的是我们之前在函数逼近中提到的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$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") #可以将看到,节点在不指定的情况下默认的是均匀样条,当然,我们可以根据散点图给#出节点的具体选择。
## 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样条基本很接近真实函数了,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$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)最后我们来讲一下怎么计算出\( 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$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)在惩罚系数为0.2的情况下,拟合还是不坏的,不是吗?至于为什么可以这样算,我们只要注意到\( \int [m^{''}(x)]dx=m^'(x_i)QR^{-1}Q^'m(x_i) \),估计的问题就与我们十分熟悉的lasso,岭回归十分相像了。本作品采用知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议进行许可。
相关文章推荐
- R语言与函数估计学习笔记(核方法与局部多项式)
- R语言与函数估计学习笔记(函数模型的参数估计)
- R语言与函数估计学习笔记(函数展开)
- [C++学习笔记]自定义函数的传值方法
- javascript学习笔记—DOM常用API、属性、方法、函数
- 黑马视频学习笔记-OC-对象方法和函数区别
- jQuery学习笔记之jQuery构建函数的7种方法
- [ADO学习笔记] Connection对象的函数与方法
- R语言与机器学习中的回归方法学习笔记
- R语言与Markov Chain Monte Carlo(MCMC)方法学习笔记(1)
- javascript学习笔记:函数与方法2
- R语言学习笔记2——常用数学函数
- JavaScript 函数的call()方法的学习笔记
- R语言与点估计学习笔记(矩估计与MLE)
- C#2005 .NET3.0高级编程学习笔记————类和结构,类的数据成员,类的函数成员(方法、属性)
- C#学习笔记------Visual C#常用函数和方法
- R语言与点估计学习笔记(刀切法与最小二乘估计)
- R语言与区间估计学习笔记
- python学习笔记之函数(方法)
- jQuery学习笔记之jQuery构建函数的7种方法