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

我对插值函数的一般求根问题感兴趣。

假设我有以下内容(x, y) data:

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 值并在图中进行注释 https://stackoverflow.com/q/37455512/4891738
• 使用拟合模型根据 Y 值预测 X 值 https://stackoverflow.com/q/43322568/4891738

建议我们简单地反转x and 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值 https://stackoverflow.com/q/52650467/4891738。该解对于线性插值来说是稳定的，但对于非线性插值来说不一定是稳定的。我现在正在寻找一个稳定的解决方案，特别是对于三次插值样条线。

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

有时，经过一段单变量线性回归y ~ x or a 单变量非线性回归y ~ f(x)我们想要反求解x对于一个目标y。这个问答就是一个例子，吸引了很多答案：求解最佳拟合多项式并绘制下拉线 https://stackoverflow.com/q/41687406/4891738，但没有一个是真正具有适应性或易于在实践中使用的。

• 接受的答案使用polyroot仅适用于简单的多项式回归；
• 使用二次公式求解析解的答案仅适用于二次多项式；
• 我的回答使用predict and uniroot一般情况下可以工作，但不方便，因为在实践中使用uniroot需要与用户交互（参见R 中的 Uniroot 解决方案 https://stackoverflow.com/q/38961221/4891738了解更多uniroot).

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

首先，让我复制中提出的线性插值的稳定解决方案我之前的回答 https://stackoverflow.com/a/52650890/4891738.

## 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" and "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
  }

It uses polyroot分段，首先找到所有根complex字段，然后只保留real分段区间上的那些。这是可行的，因为三次插值样条只是许多分段三次多项式。我的回答在如何在R中保存和加载样条插值函数？ https://stackoverflow.com/q/51540182/4891738已经展示了如何获得分段多项式系数，因此使用polyroot很简单。

使用问题中的示例数据，两者RootSpline1 and 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