获取给定 y 值的 x 值：线性/非线性插值函数的一般求根

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 值并在图中对其进行注释
• 使用拟合模型从 Y 值预测 X 值

建议我们简单地反转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

首先，根数不对。我们从图中（左侧）看到两个根，但代码只返回一个。其次，它不是一个正确的根，因为

f1(fake_root)
#[1] 2.906103

不是 2.85。

我在如何在 R 中的 approxfun() 之后从 y 值输入估计 x 值中对这个一般问题进行了第一次尝试。该解对于线性插值是稳定的，但对于非线性插值不一定稳定。我现在正在寻找一个稳定的解决方案，特别是三次插值样条。

---

解决方案如何在实践中发挥作用？

有时在单变量线性回归y ~ x或单变量非线性回归之后，y ~ f(x)我们想要对x目标进行反向求解y。这个问答是一个例子，吸引了很多答案：求解最佳拟合多项式和绘制下拉线，但没有一个是真正自适应的或在实践中易于使用。

• 接受的答案 usingpolyroot仅适用于简单的多项式回归；
• 使用二次公式作为解析解的答案仅适用于二次多项式；
• 我的答案使用predict和uniroot一般有效，但并不方便，因为在实践中使用uniroot需要与用户交互（有关更多信息，请参阅R 中的 Uniroot 解决方案uniroot）。

如果有一个自适应且易于使用的解决方案，那就太好了。

Answer 1

首先，让我复制在我之前的答案中提出的线性插值的稳定解决方案。

## 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分段，首先在复数域上找到所有的根，然后在分段区间上只保留实数。这是有效的，因为三次插值样条只是一些分段三次多项式。我对如何在 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