是否有非结三次样条函数?

Is there a function for not-a-knot cubic splines?

我需要在我的 R 脚本中使用“非结”三次样条进行插值。

虽然有一些样条 R 包,但其中 none 似乎考虑了“非结”类型,即使据说它是一种相当“流行”的三次样条类型, 它在 Matlab 中可用。

我担心“非结”三次样条还有另一个名称。它是一个三次样条,其中两个额外条件是关于第二个和最后一个节点中的三阶导数连续性(而不是像自然三次样条或其他选择那样将一阶导数固定在端点节点处)。

稍微挖掘一下,似乎是在pracma::interp1

中实现的

来自?pracma::interp1

Interpolation to find ‘yi’, the values of the underlying function at the points in the vector ‘xi’. ... Method `spline' uses the spline approach by Moler et al., and is identical with the Matlab option of the same name, but slightly different from R's spline function.

这里没有提到“not-a-knot”,但我是通过使用 sos::findFn("not-a-knot") 到达那里的,这将我带到了 gsignal::pulstran()。该函数的“样条”方法被描述为“使用 not-a-knot 结束条件的样条插值”; method= 参数直接传递给 pracma::interp1.

这是一个比较的例子:

  • 两个 base-R 样条实现
  • pracma::interp1
  • pracma:::.ml.spline(不做任何argument-checking,排序,去重...但是允许外推。我不知道为什么pracma::interp1做这个限制。软件包已开发 on r-forge 但我没有看到邮件列表等,您可以联系维护者或自己破解 interp1 ...)
  • Octave(改编自 :见下面的代码。有一个 RcppOctave 包,但它已存档,我无法通过 remotes::install_version("RcppOctave", version = "0.18-1") 轻松获得要安装的包 -配置漂移太多...)
library(splines)
library(pracma)
x <- 1:6
y <- c(16, 18, 21, 17, 15, 12)
xi0 <- seq(1, 6, by = 0.1)
xi1 <- seq(1, 7, by = 0.1)
ys1 <- interp1(x, y, xi0, method = "spline")
ys2 <- predict(interpSpline(x, y), xi1)$y
ys3 <- spline(x, y, xout  = xi1)$y
ys4 <- pracma:::.ml.spline(x, y, xi1)
ys5 <- scan("octave_out.mat", comment = "#") ## see below

png("spline2.png")
par(las=1, bty = "l", cex = 1.5, lwd = 2)
plot(x,y, xlim = range(xi1), ylim = range(ys4),
     pch = 16, col = "gray")
cvec <- c(palette()[c(1,2,4,6)],
  ## make last color semi-transparent so we can see the overlay
          adjustcolor(palette()[5], alpha.f = 0.4))
matlines(xi1, cbind(ys2,ys3,ys4, ys5) , col = cvec[1:4], lwd = 2)
lines(xi0, ys1, lwd = 3, col = cvec[5])
legend("bottomleft",
       lty = c(1:4, 1),
       col = cvec,
       legend = c("interpSpline", "spline", ".ml.spine", "octave", "interp1"))
dev.off()

  • Octave 和 .ml.spline 输出一致
  • interp1 在数据范围内一致(无外推)
  • 两个 R 样条给出了不同但合理的答案

八度代码(运行在一个单独的shell)

x = 1:6;
y = [16, 18, 21, 17, 15, 12];
xi1 = linspace(1, 7, 61);
pp = interp1(x, y, 'spline', 'pp');
yy = ppval(pp, xi1);
save octave_out.mat yy