R语言与函数估计学习笔记(样条方法)
样条估计
如果函数在不同地方有不同的非线性度,或者有多个极值点,那么用多项式特别是低阶多项式来完成拟合是非常不合适的。一种解决办法是我们之前提到的近邻多项式(或者称局部多项式),另一种就是样条——用分段的低阶多项式逼近函数。
关于样条,常用的有两类,一类是多项式样条,另一类是光滑样条。
多项式样条
多项式样条的样条基有很多,最为著名的是我们之前在函数逼近中提到的truncated power basis与B-spline basis。我们这里十分简要的介绍一下B样条,B样条基下的函数逼近可以写为:
其中
上式中否则取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)满足如下的惩罚最小二乘:
在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+λQ′Q使得它可以进行Cholesky分解,有因子L,DStep 3: 解方程:(R+λQ′Q)γ=Q′yStep 4: 得到估值m=y−αQγ.上面的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,岭回归十分相像了。
数据分析咨询请扫描二维码
若不方便扫码,搜微信号:CDAshujufenxi
在数字化浪潮中,数据驱动决策已成为企业发展的核心竞争力,数据分析人才的需求持续飙升。世界经济论坛发布的《未来就业报告》, ...
2025-03-28你有没有遇到过这样的情况?流量进来了,转化率却不高,辛辛苦苦拉来的用户,最后大部分都悄无声息地离开了,这时候漏斗分析就非 ...
2025-03-27TensorFlow Datasets(TFDS)是一个用于下载、管理和预处理机器学习数据集的库。它提供了易于使用的API,允许用户从现有集合中 ...
2025-03-26"不谋全局者,不足谋一域。"在数据驱动的商业时代,战略级数据分析能力已成为职场核心竞争力。《CDA二级教材:商业策略数据分析 ...
2025-03-26当你在某宝刷到【猜你喜欢】时,当抖音精准推来你的梦中情猫时,当美团外卖弹窗刚好是你想吃的火锅店…… 恭喜你,你正在被用户 ...
2025-03-26当面试官问起随机森林时,他到底在考察什么? ""请解释随机森林的原理""——这是数据分析岗位面试中的经典问题。但你可能不知道 ...
2025-03-25在数字化浪潮席卷的当下,数据俨然成为企业的命脉,贯穿于业务运作的各个环节。从线上到线下,从平台的交易数据,到门店的运营 ...
2025-03-25在互联网和移动应用领域,DAU(日活跃用户数)是一个耳熟能详的指标。无论是产品经理、运营,还是数据分析师,DAU都是衡量产品 ...
2025-03-24ABtest做的好,产品优化效果差不了!可见ABtest在评估优化策略的效果方面地位还是很高的,那么如何在业务中应用ABtest? 结合企业 ...
2025-03-21在企业数据分析中,指标体系是至关重要的工具。不仅帮助企业统一数据标准、提升数据质量,还能为业务决策提供有力支持。本文将围 ...
2025-03-20解锁数据分析师高薪密码,CDA 脱产就业班助你逆袭! 在数字化浪潮中,数据驱动决策已成为企业发展的核心竞争力,数据分析人才的 ...
2025-03-19在 MySQL 数据库中,查询一张表但是不包含某个字段可以通过以下两种方法实现:使用 SELECT 子句以明确指定想要的字段,或者使 ...
2025-03-17在当今数字化时代,数据成为企业发展的关键驱动力,而用户画像作为数据分析的重要成果,改变了企业理解用户、开展业务的方式。无 ...
2025-03-172025年是智能体(AI Agent)的元年,大模型和智能体的发展比较迅猛。感觉年初的deepseek刚火没多久,这几天Manus又成为媒体头条 ...
2025-03-14以下的文章内容来源于柯家媛老师的专栏,如果您想阅读专栏《小白必备的数据思维课》,点击下方链接 https://edu.cda.cn/goods/sh ...
2025-03-13以下的文章内容来源于刘静老师的专栏,如果您想阅读专栏《10大业务分析模型突破业务瓶颈》,点击下方链接 https://edu.cda.cn/go ...
2025-03-12以下的文章内容来源于柯家媛老师的专栏,如果您想阅读专栏《小白必备的数据思维课》,点击下方链接 https://edu.cda.cn/goods/sh ...
2025-03-11随着数字化转型的加速,企业积累了海量数据,如何从这些数据中挖掘有价值的信息,成为企业提升竞争力的关键。CDA认证考试体系应 ...
2025-03-10推荐学习书籍 《CDA一级教材》在线电子版正式上线CDA网校,为你提供系统、实用、前沿的学习资源,助你轻松迈入数据分析的大门! ...
2025-03-07在数据驱动决策的时代,掌握多样的数据分析方法,就如同拥有了开启宝藏的多把钥匙,能帮助我们从海量数据中挖掘出关键信息,本 ...
2025-03-06