R smooth.spline 拟合平滑样条曲线

说明

将三次平滑样条拟合到提供的数据。

用法

smooth.spline(x, y = NULL, w = NULL, df, spar = NULL, lambda = NULL, cv = FALSE,
              all.knots = FALSE, nknots = .nknots.smspl,
              keep.data = TRUE, df.offset = 0, penalty = 1,
              control.spar = list(), tol = 1e-6 * IQR(x), keep.stuff = FALSE)

.nknots.smspl(n)

参数

x	给出预测变量值的向量，或指定 x 和 y 的列表或两列矩阵。
y	回应。如果 y 丢失或 NULL ，则假定响应由 x 指定，其中 x 为索引向量。
w	与 x 长度相同的可选权重向量；默认全部为1。
df	所需的等效自由度数(平滑矩阵的迹)。必须在 (1,n_x] 、 n_x 中唯一 x 值的数量，请参见下文。
spar	平滑参数，通常(但不一定)位于 (0,1] 中。当指定 spar 时，拟合(惩罚对数似然)准则中的二阶导数平方积分的系数 \lambda 是 spar 的单调函数，请参阅下面的详细信息。或者，可以指定 lambda 来代替无标度 spar = s 。
lambda	如果需要，可以指定内部(取决于设计)平滑参数 \lambda 而不是 spar 。这对于交叉验证或引导程序等重采样算法可能是理想的。
cv	当 FALSE 时，普通留一法 ( TRUE ) 或 ‘generalized’ 交叉验证 (GCV) ；仅当spar和df均未指定时才用于平滑参数计算；然而，它用于确定结果中的cv.crit。将其设置为 NA 以获得加速会跳过杠杆和任何分数的评估。
all.knots	如果 TRUE ，则 x 中的所有不同点都用作结。如果 FALSE (默认)，则使用 x[] 的子集，特别是 x[j]，其中 nknots 索引在 1:n 中均匀分布，另请参阅下一个参数 nknots 。 或者，使用严格递增的numeric向量，指定要使用的“all the knots”；必须重新缩放为 [0, 1]，使其对应于返回的 ans $ fit$knots 序列，而不重复边界结。
nknots	整数或 function 给出 all.knots = FALSE 时要使用的结数。如果是函数(默认情况下)，则结数为 nknots(nx) 。默认情况下使用 .nknots.smspl() ，对于 n_x > 49 这小于 n_x ，即唯一 x 值的数量，请参阅注释。
keep.data	逻辑指定输入数据是否应保留在结果中。如果TRUE(默认)，则可从结果中获得拟合值和残差。
df.offset	允许 GCV 准则中的自由度增加df.offset。
penalty	GCV 准则中自由度的惩罚系数。
control.spar	计算平滑参数 spar 时控制根查找的命名组件的可选列表，即缺少或 NULL ，请参见下文。 请注意，这部分是实验性的，可能会随着一般晶石计算的改进而改变！ 低的： 下界为spar;默认为 -1.5(用于隐式默认为 0R1.4 之前的版本)。 高的： spar 的上限；默认为+1.5。 托尔： 使用的绝对精度(公差)；默认为 1e-4(以前为 1e-3)。 每股收益： 使用的相对精度；默认为 2e-8(以前为 0.00244)。 痕迹： 逻辑指示是否应该跟踪迭代。 马克西特： 给出最大迭代次数的整数；默认为 500。 请注意，仅在区间 [low, high] 中搜索 spar 。
tol	same-ness 的容差或 x 值的唯一性。这些值被分箱到大小为 tol 的箱中，并且落入同一箱中的值被视为相同。必须严格为正(且有限)。
keep.stuff	实验性的 logical 指示结果是否应保留内部计算的额外内容。应该允许重建 X 矩阵等。
n	对于.nknots.smspl；通常是唯一 x 值(又名 n_x )的数量。

细节

x 和 y 都不允许包含缺失值或无限值。

x 向量应至少包含四个不同的值。这里的‘Distinct’由tol控制：被视为相同的值被它们的第一个值替换，并且相应的y和w被相应地池化。

除非指定 lambda 而不是 spar ，否则使用的计算 \lambda (作为 s=spar 的函数)为 \lambda = r \cdot 256^{3 s - 1} ，其中 r = tr(X' W X) / tr(\Sigma) 、 \Sigma 是由 \Sigma_{ij} = \int B_i''(t) B_j''(t) dt 给出的矩阵， X 由 X_{ij} = B_j(x_i) 给出，W 是权重的对角矩阵(缩放后其轨迹为 n，即原始观测数)，B_k(.) 是 k -th B-spline 。

请注意，使用这些定义 f_i = f(x_i) 和 B-spline 基础表示 f = X c (即 c 是样条系数的向量)，惩罚对数似然为 L = (y - f)' W (y - f) + \lambda c' \Sigma c ，因此 c 是(岭回归)(X' W X + \lambda \Sigma) c = X' W y 的解决方案。

如果 spar 和 lambda 缺失或 NULL ，则使用 df 的值来确定平滑程度。如果 df 也丢失，则使用留一交叉验证(普通或 ‘generalized’ 由 cv 确定)来确定 \lambda 。

请注意，根据上述关系，spar是s = s0 + 0.0601 \cdot \log\lambda，这是故意的不同的来自 S-PLUS 的实现smooth.spline(在哪里spar正比于\lambda)。在R的(\log \lambda)规模，改变更有意义spar线性地。

但请注意，目前对于小于 -1 或 -2 的 spar 值，结果可能变得非常不可靠。对于大于 2 左右的值也可能会发生同样的情况。不要考虑将 spar 或控件 low 和 high 设置在这样的安全范围之外，除非您知道自己在做什么！同样，指定 lambda 而不是 spar 也很微妙，特别是 lambda 的 “safe” 值的范围不是比例不变的，因此完全依赖于数据。

当 x 中存在重复点时，‘generalized’ 交叉验证方法 GCV 将正常工作。然而，留一交叉验证对于重复点的含义并不明确，并且内部代码使用涉及省略重复点组的近似值。在这种情况下最好避免cv = TRUE。

值

带有组件的 "smooth.spline" 类的对象

x	不同的 x 值按升序排列，请参阅上面的“详细信息”。
y	对应于 x 的拟合值。
w	在 x 的唯一值处使用的权重。
yin	用于唯一 y 值的 y 值。
tol	tol 参数(其默认值取决于 x )。
data	仅当keep.data = TRUE：本身是一个list，其组件x、y和w具有相同的长度。这些是原始的 (x_i,y_i,w_i), i = 1, \dots, n 值，其中 data$x 可能具有重复值，因此比上述 x 组件长；查看具体信息。
n	一个整数； (原始)样本量。
lev	(当 cv 不是 NA 时)利用平滑矩阵的对角线值。
cv	使用的 cv 参数；即 FALSE 、 TRUE 或 NA 。
cv.crit	交叉验证分数，‘generalized’ 或 true，具体取决于 cv 。 CV 分数通常称为 “PRESS”(并标记为 print() )，表示“预测平方和”。请注意，这与拟合期间最小化的(CV 或 GCV)分数不同(并在 crit 中返回)，例如，在 nx < n 的情况下(其中 nx =n_x 是唯一的 x 值)。
pen.crit	惩罚标准，一个非负数；只是(加权)残差平方和 (RSS)， sum(.$w * residuals(.)^2) 。
crit	底层标准值最小化.Fortran常规 'sslvrg’。什么时候df已指定，标准为3 + (tr(S_\lambda) - df)^2, 其中3 +是出于数字(和历史)原因吗？
df	使用等效自由度。请注意，(目前)当真实的 df 介于 1 和 2 之间时，该值可能会变得相当不精确。
spar	计算或给出的 spar 的值，除非在此处设置为 NA 时将其指定为 c(lambda = *) 。
ratio	(当上面的 spar 不是 NA 时)，值 r ，两个矩阵迹线的比率。
lambda	\lambda 对应 spar 的值，请参见上面的详细信息。
iparms	命名的整数(3)向量，其中..$ipars["iter"]给出了使用的spar计算迭代次数。
auxMat	实验性的；当 keep.stuff 为 true 时，包含部分内部计算的 “flat” 数值向量。
fit	供 predict.smooth.spline 使用的列表，包含组件 结： 结序列(包括重复的边界结)，缩放为 [0, 1] (通过 min 和 range )。 NK: 系数数或 ‘proper’ 结数加 2。 系数： 使用的样条基础的系数。 最小值，范围： 给出 x 相应数量的数字。
call	匹配的调用。

method(class = "smooth.spline") 显示了基于上面lev 向量的hatvalues() 方法。

注意

唯一 x 值的数量 \mathbf{nx} = n_x 由 tol 参数确定，相当于

    nx <- length(x) - sum(duplicated( round((x - mean(x)) / tol) ))
  

默认all.knots = FALSE和nknots = .nknots.smspl，只需要使用O({n_x}^{0.2})结而不是n_x为了n_x > 49。这降低了速度和内存需求，但不再大幅降低，因为R版本 1.5.1 仅适用于O(n_k) + O(n)其中n_k是节数。

在这种情况下，并非所有唯一的 x 值都用作结，结果是回归样条而不是严格意义上的平滑样条，但非常接近，除非使用较小的平滑参数(或较大的 df )。

例子

require(graphics)
plot(dist ~ speed, data = cars, main = "data(cars)  &  smoothing splines")
cars.spl <- with(cars, smooth.spline(speed, dist))
cars.spl
## This example has duplicate points, so avoid cv = TRUE

lines(cars.spl, col = "blue")
ss10 <- smooth.spline(cars[,"speed"], cars[,"dist"], df = 10)
lines(ss10, lty = 2, col = "red")
legend(5,120,c(paste("default [C.V.] => df =",round(cars.spl$df,1)),
               "s( * , df = 10)"), col = c("blue","red"), lty = 1:2,
       bg = 'bisque')


## Residual (Tukey Anscombe) plot:
plot(residuals(cars.spl) ~ fitted(cars.spl))
abline(h = 0, col = "gray")

## consistency check:
stopifnot(all.equal(cars$dist,
                    fitted(cars.spl) + residuals(cars.spl)))
## The chosen inner knots in original x-scale :
with(cars.spl$fit, min + range * knot[-c(1:3, nk+1 +1:3)]) # == unique(cars$speed)

## Visualize the behavior of  .nknots.smspl()
nKnots <- Vectorize(.nknots.smspl) ; c.. <- adjustcolor("gray20",.5)
curve(nKnots, 1, 250, n=250)
abline(0,1, lty=2, col=c..); text(90,90,"y = x", col=c.., adj=-.25)
abline(h=100,lty=2); abline(v=200, lty=2)

n <- c(1:799, seq(800, 3490, by=10), seq(3500, 10000, by = 50))
plot(n, nKnots(n), type="l", main = "Vectorize(.nknots.smspl) (n)")
abline(0,1, lty=2, col=c..); text(180,180,"y = x", col=c..)
n0 <- c(50, 200, 800, 3200); c0 <- adjustcolor("blue3", .5)
lines(n0, nKnots(n0), type="h", col=c0)
axis(1, at=n0, line=-2, col.ticks=c0, col=NA, col.axis=c0)
axis(4, at=.nknots.smspl(10000), line=-.5, col=c..,col.axis=c.., las=1)

##-- artificial example
y18 <- c(1:3, 5, 4, 7:3, 2*(2:5), rep(10, 4))
xx  <- seq(1, length(y18), length.out = 201)
(s2   <- smooth.spline(y18)) # GCV
(s02  <- smooth.spline(y18, spar = 0.2))
(s02. <- smooth.spline(y18, spar = 0.2, cv = NA))
plot(y18, main = deparse(s2$call), col.main = 2)
lines(s2, col = "gray"); lines(predict(s2, xx), col = 2)
lines(predict(s02, xx), col = 3); mtext(deparse(s02$call), col = 3)

## Specifying 'lambda' instead of usual spar :
(s2. <- smooth.spline(y18, lambda = s2$lambda, tol = s2$tol))



## The following shows the problematic behavior of 'spar' searching:
(s2  <- smooth.spline(y18, control =
                      list(trace = TRUE, tol = 1e-6, low = -1.5)))
(s2m <- smooth.spline(y18, cv = TRUE, control =
                      list(trace = TRUE, tol = 1e-6, low = -1.5)))
## both above do quite similarly (Df = 8.5 +- 0.2)

作者

R implementation by B. D. Ripley and Martin Maechler (spar/lambda, etc).

来源

该函数基于 T. Hastie 和 R. Tibshirani 编写的 GAMFIT Fortran 程序中的代码(最初取自 http://lib.stat.cmu.edu/general/gamfit )，该程序使用 Finbarr O'Sullivan 的样条代码。其设计与 Chambers & Hastie (1992) 的 smooth.spline 函数类似。

参考

Chambers, J. M. and Hastie, T. J. (1992) Statistical Models in S, Wadsworth & Brooks/Cole.

Green, P. J. and Silverman, B. W. (1994) Nonparametric Regression and Generalized Linear Models: A Roughness Penalty Approach. Chapman and Hall.

Hastie, T. J. and Tibshirani, R. J. (1990) Generalized Additive Models. Chapman and Hall.

也可以看看

predict.smooth.spline 用于评估样条曲线及其导数。