получить значение x по заданному значению y: общий корень для линейной / нелинейной интерполяционной функции

Question

Меня интересует общая проблема поиска корней для функции интерполяции.

Предположим, у меня есть следующие данные (x, y):

set.seed(0)
x <- 1:10 + runif(10, -0.1, 0.1)
y <- rnorm(10, 3, 1)

, а также линейная интерполяция иинтерполяция кубического сплайна:

f1 <- approxfun(x, y)
f3 <- splinefun(x, y, method = "fmm")

Как найти значения x, в которых эти функции интерполяции пересекают горизонтальную линию y = y0?Ниже приведена графическая иллюстрация с y0 = 2.85.

par(mfrow = c(1, 2))
curve(f1, from = x[1], to = x[10]); abline(h = 2.85, lty = 2)
curve(f3, from = x[1], to = x[10]); abline(h = 2.85, lty = 2)

Мне известны несколько предыдущих тем по этой теме, например

• прогнозирование значений x по простой подгонке и аннотирование их на графике
• прогнозирование значения X по значению Y с помощью подгоночной модели

Предлагается просто поменять местами x и y, выполнить интерполяцию для (y, x) и вычислить интерполированное значение в y = y0.

Однако это фальшивкаидея.Пусть y = f(x) будет интерполяционной функцией для (x, y), эта идея верна только тогда, когда f(x) является монотонной функцией x, так что f является обратимым.В противном случае x не является функцией y и интерполяция (y, x) не имеет смысла.

Принимая линейную интерполяцию с моими примерами, эта ложная идея дает

fake_root <- approx(y, x, 2.85)[[2]]
# [1] 6.565559

Firstколичество корней неверно.Мы видим два корня из рисунка (слева), но код возвращает только один.Во-вторых, это не правильный корень, так как

f1(fake_root)
#[1] 2.906103

не равен 2,85.

Моя первая попытка решить эту общую проблему: Как оценить значение x по yввод значения после прибл. () в R .Решение оказывается устойчивым для линейной интерполяции, но не обязательно устойчивым для нелинейной интерполяции.Сейчас я ищу стабильное решение, особенно для кубического сплайн-сплайна.

---

Как решение может быть полезным на практике?

Иногда после одномерного линейная регрессия y ~ x или одномерный нелинейная регрессия y ~ f(x), для которой мы хотим выполнить обратное решение x для цели y.Эти вопросы и ответы являются примером и привлекли много ответов: Решите наиболее подходящие полиномиальные линии и выпадающие линии графика , но ни одна из них не является действительно адаптивной или простой в использовании на практике.

• Принятый ответ с использованием polyroot работает только для простой полиномиальной регрессии;
• Ответы с использованием квадратичной формулы для аналитического решения работают только для квадратичного полинома;
• Мой ответ с использованием predict иuniroot работает в целом, но не удобно, так как на практике использование uniroot требует взаимодействия с пользователями (подробнее см. Решение Uniroot в R для uniroot).

Было бы очень хорошо, если бы существовало адаптивное и простое в использовании решение.

Answer 1

С учетом точек данных и функции сплайна, как указано выше, просто примените findzeros() из пакета pracma .

library(pracma)
xs <- findzeros(function(x) f3(x) - 2.85,min(x), max(x))

xs  # [1] 3.924513 6.435812 9.207169 9.886618
points(xs, f3(xs))

Answer 2

Прежде всего, позвольте мне скопировать в стабильное решение для линейной интерполяции, предложенное в мой предыдущий ответ .

## given (x, y) data, find x where the linear interpolation crosses y = y0
## the default value y0 = 0 implies root finding
## since linear interpolation is just a linear spline interpolation
## the function is named RootSpline1
RootSpline1 <- function (x, y, y0 = 0, verbose = TRUE) {
  if (is.unsorted(x)) {
     ind <- order(x)
     x <- x[ind]; y <- y[ind]
     }
  z <- y - y0
  ## which piecewise linear segment crosses zero?
  k <- which(z[-1] * z[-length(z)] <= 0)
  ## analytical root finding
  xr <- x[k] - z[k] * (x[k + 1] - x[k]) / (z[k + 1] - z[k])
  ## make a plot?
  if (verbose) {
    plot(x, y, "l"); abline(h = y0, lty = 2)
    points(xr, rep.int(y0, length(xr)))
    }
  ## return roots
  xr
  }

Для сплайнов кубической интерполяции, возвращаемых stats::splinefun методами "fmm", "natrual", "periodic" и "hyman", следующая функция обеспечивает устойчивое численное решение.

RootSpline3 <- function (f, y0 = 0, verbose = TRUE) {
  ## extract piecewise construction info
  info <- environment(f)$z
  n_pieces <- info$n - 1L
  x <- info$x; y <- info$y
  b <- info$b; c <- info$c; d <- info$d
  ## list of roots on each piece
  xr <- vector("list", n_pieces)
  ## loop through pieces
  i <- 1L
  while (i <= n_pieces) {
    ## complex roots
    croots <- polyroot(c(y[i] - y0, b[i], c[i], d[i]))
    ## real roots (be careful when testing 0 for floating point numbers)
    rroots <- Re(croots)[round(Im(croots), 10) == 0]
    ## the parametrization is for (x - x[i]), so need to shift the roots
    rroots <- rroots + x[i]
    ## real roots in (x[i], x[i + 1])
    xr[[i]] <- rroots[(rroots >= x[i]) & (rroots <= x[i + 1])]
    ## next piece
    i <- i + 1L
    }
  ## collapse list to atomic vector
  xr <- unlist(xr)
  ## make a plot?
  if (verbose) {
    curve(f, from = x[1], to = x[n_pieces + 1], xlab = "x", ylab = "f(x)")
    abline(h = y0, lty = 2)
    points(xr, rep.int(y0, length(xr)))
    }
  ## return roots
  xr
  }

Она использует polyroot кусочно, сначала находя все корни в complex *Поле 1016 *, затем на кусочном интервале сохраняются только действительные единицы.Это работает, потому что кубический интерполяционный сплайн - это просто количество кусочно-кубических полиномов.Мой ответ на Как сохранить и загрузить сплайн-интерполяционные функции в R? показал, как получить кусочно-полиномиальные коэффициенты, поэтому использование polyroot не вызывает затруднений.

Использование примеров данных в вопросеи RootSpline1 и RootSpline3 правильно идентифицируют все корни.

par(mfrow = c(1, 2))
RootSpline1(x, y, 2.85)
#[1] 3.495375 6.606465
RootSpline3(f3, 2.85)
#[1] 3.924512 6.435812 9.207171 9.886640