интегрировать квадратичные b-сплайны в R

Я работаю с функцией, которая зависит от квадратичной интерполяции B-сплайна, предварительно оцененной функцией cobs в том же пакете R. Расчетные узлы и соответствующие коэффициенты приведены в коде. Далее мне нужен интеграл этой функции от 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 integrate spline bspline
person Akkariz    schedule 04.02.2018    source источник

Ответы (1)


arrow_upward
2
arrow_downward

Я думаю, что алгоритм, используемый функцией «интегрировать», не подходит для этих условий. Например, если вы измените нижние пределы, все будет работать так, как ожидалось:

> 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)
person Esteban PS    schedule 04.02.2018