在 R 中获取 smooth.splines 的公式答案

【问题标题】：Getting Formulas for smooth.splines in R在 R 中获取 smooth.splines 的公式
【发布时间】：2021-12-08 21:53:19
【问题描述】：

我使用了 stats::smooth.spline 函数来拟合 60 (x, y) 对的数据集，现在我需要访问可以为我生成预测的公式，同时我可以访问模型。不幸的是，文档 here 并没有帮助我解决这个问题（而且看起来也有点过时了）。

我的理解是smooth.spline 提供了一个函数g: Reals -> Reals，其中g 是满足g = argmin(g) spar * SUM(MSE) + lambda * Integral [g''(x)]^2 dx 的三次多项式。如果是这种情况，我希望看到 g(x) = x^3 + 5x^2 + 10x + 15 的输出（尽管生成的对象似乎只包含 spar 和 lambda）。

smooth.spline 也有可能沿着数据的不同部分形成多个样条，在这种情况下，我想要的输出类似于：

在 x=0 到 x=2 上，预测 y = 3x^3 + 2x^2 + x
在 x=2 到 x=3.5 上，预测 y = x^3 + 2x^2 + 3x + 12
......

下面是一些示例代码：

dataset <- data.frame(x = c(1,2,3,4,5,6,7), y = c(10,8,4,6,8,11,15))
spline_result <- smooth.spline(dataset$x, y = dataset$y)
# Plotting the spline looks like a parabola centered near x=3.5 y=5
# so I'd expect something like g(x) = k *(x - 3.5) ^ 2 + 5
# where k is some constant

附加信息：我正在运行 R 版本 4.1.2 2021-11-01 Bird Hippie with packageVersion("stats") 4.1.2

感谢您的帮助！

【问题讨论】：

smooth.spline 在每个唯一的 x 值之间使用不同的三次多项式，并且在 x 值的范围之外是线性的。要找到系数，您可以对每个区间进行 4 次预测并使用多项式回归来拟合值，但这不太可能是一个有用的练习。

标签： r spline cubic-spline

【解决方案1】：

smooth.spline 函数在每个x 值处生成一个带有结的“自然样条曲线”。这意味着它是x 范围之外的线性函数，并且是每个值之间的三次多项式。因此，在您的情况下，您将有 2 个线性方程和 6 个三次方程。这通常不是使用此类函数的有用方法，因为系数通常几乎相互抵消，因此当您尝试评估它时可能会出现很多舍入误差。

如果您只是对评估不在原始x 向量中的各个点的样条曲线感兴趣，请使用predict() 函数，例如

dataset <- data.frame(x = c(1,2,3,4,5,6,7), y = c(10,8,4,6,8,11,15))
spline_result <- smooth.spline(dataset$x, y = dataset$y)

newx <- seq(0, 10, len = 100)
newvals <- predict(spline_result, x = newx)
plot(newvals, type = 'l')
points(dataset)

^{由reprex package (v2.0.1) 于 2021-12-09 创建}

predict() 函数通过避免多项式的幂基础来避免舍入误差。

如果您真的想要多项式系数，获得它们的一种方法是对预测使用多项式回归。例如，要查找 3 到 4 之间的段的系数，您可以使用

lm(y ~ poly(x, degree = 3, raw = TRUE), data = predict(spline_result, x = seq(3, 4, len = 10)))
#> 
#> Call:
#> lm(formula = y ~ poly(x, degree = 3, raw = TRUE), data = predict(spline_result, 
#>     x = seq(3, 4, len = 10)))
#> 
#> Coefficients:
#>                      (Intercept)  poly(x, degree = 3, raw = TRUE)1  
#>                          26.3378                          -14.9943  
#> poly(x, degree = 3, raw = TRUE)2  poly(x, degree = 3, raw = TRUE)3  
#>                           3.2874                           -0.2059

^{由reprex package (v2.0.1) 于 2021-12-09 创建}

这给出多项式 26.3378 -14.9943 x + 3.2874 x^2 -0.2059 x^3。

【讨论】：