在给定 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")

如何找到这些插值函数穿过水平线 y = y0 的 x 值？下面是带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简单反转，对(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

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

f1(fake_root)
#[1] 2.906103

不是 2.85。

我在第一次尝试解决这个一般性问题。该解决方案对于线性插值是稳定的，但对于非线性插值则不一定稳定。我现在正在寻找一个稳定的解决方案，专门针对三次插值样条。

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

有时在单变量线性回归y ~ x或单变量非线性回归y ~ f(x)之后，我们想要为目标 y 反向求解 x。这个问答是一个例子，吸引了很多回答：，但是none是真正适应性强或者在实践中好用的

接受的答案使用 polyroot 仅适用于简单的多项式回归；
使用二次公式求解解析解的答案仅适用于二次多项式；
我使用 predict 和 uniroot 的答案一般有效，但不方便，因为在实践中使用 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 分段，首先在 complex 字段上找到所有根，然后在分段上只保留 real间隔。这是有效的，因为三次插值样条只是一些分段三次多项式。我在的回答已经展示了如何获得分段多项式系数，因此使用 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

Answer 2

给定数据点和样条函数，只需应用 pracma 包中的 findzeros()。

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))

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

get x-value given y-value: general root finding for linear / non-linear interpolation function

r

regression

interpolation

spline

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