1

我正在使用一个函数，该函数依赖于由cobs同一R包中的函数预先估计的二次 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函数上集成，并相应地转换结果。但是，我很想知道为什么会发生这种奇怪的行为。

4

1 回答 1

2

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

> 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

这在数值积分方法中很常见，您必须根据具体情况进行选择。我正在使用梯形规则在同一区域上进行集成并且运行良好的原始代码

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)
于 2018-02-04T10:50:03.477 回答