在 R 中积分二次 b 样条答案

【问题标题】：integrate quadratic b-splines in R在 R 中积分二次 b 样条
【发布时间】：2018-02-04 10:16:30
【问题描述】：

我正在使用一个函数，该函数依赖于由同一 R 包中的 cobs 函数预先估计的二次 B 样条插值。估计的结和相应的系数在代码中给出。更进一步，我需要这个函数从 0 到某个值的积分，例如 0.6 或 0.7。由于我的函数是严格正数，如果积分上限增加，积分值应该增加。但是，对于某些值，情况并非如此，如使用 0.6 和 0.7 时所示

library(cobs)
b <- 0.6724027
xi1 <- 0.002541667
xi2 <- 2.509625
knots <- c(5.000010e-06, 8.700000e-05, 3.420000e-04, 1.344000e-03, 5.292000e-03, 2.082900e-02, 8.198800e-02, 3.227180e-01, 1.270272e+00, 5.000005e+00)
coef <- c(2.509493, 2.508141, 2.466733, 2.378368, 2.239769, 2.063977, 1.874705, 1.601780, 1.288163, 1.262683, 1.432729)
fn <- function(x) {
  z <- (2 - b) * (cobs:::.splValue(2, knots, coef, x, 0) - 2 * x * xi1) / xi2 - b
  return (z)
}

x <- seq(0, 0.7, 0.0001)
plot(x, fn(x), type = 'l')
integrate(f = fn, 0, 0.6)
# 0.1049019 with absolute error < 1.2e-15
integrate(f = fn, 0, 0.7)
# 0.09714124 with absolute error < 1.1e-15

我知道我可以直接在cobs:::.splValue 函数上集成，并相应地转换结果。但是，我很想知道为什么会发生这种奇怪的行为。

【问题讨论】：

标签： r spline integrate bspline

【解决方案1】：

我认为函数“integrate”使用的算法在这些条件下表现不佳。例如，如果您修改下限，它会按预期工作：

> integrate(f = fn, 0.1, 0.6)
0.06794357 with absolute error < 7.5e-16

> integrate(f = fn, 0.1, 0.7)
0.07432096 with absolute error < 8.3e-16

这在数值积分方法中很常见，您必须根据具体情况进行选择。我正在使用梯形规则在同一区域上进行积分并且效果很好original code

composite.trapezoid <- function(f, a, b, n) {
  if (is.function(f) == FALSE) {
    stop('f must be a function with one parameter (variable)')
  }

  h <- (b - a) / n

  j <- 1(:n - 1)

  xj <- a + j * h

  approx <- (h / 2) * (f(a) + 2 * sum(f(xj)) + f(b))

  return(approx)
}

> composite.trapezoid(f = fn, 0, 0.6, 10000)
[1] 0.1079356
> composite.trapezoid(f = fn, 0, 0.7, 10000)
[1] 0.1143195

如果我们分析积分接近0.65区域的行为，我们可以看到第一种方法存在问题（不平滑）：

tst = sapply(seq(0.5, 0.8, length.out = 100), function(upper) {
  integrate(f = fn, 0, upper)[[1]]

})
plot(seq(0.5, 0.8, length.out = 100), tst)

并且梯形规则表现更好：

tst2 = sapply(seq(0.5, 0.8, length.out = 100), function(upper) {
  composite.trapezoid(f = fn, 0, upper, 10000)[[1]]

})
plot(seq(0.5, 0.8, length.out = 100), tst2)

【讨论】：