在给定y值的情况下获得x值:线性/非线性插值函数的一般根查找

问题描述

我对插值函数的一般根查找问题感兴趣.

I am interested in a general root finding problem for an interpolation function.

假设我有以下(x, y)数据:

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

以及线性插值和三次样条插值:

as well as a linear interpolation and a cubic spline interpolation:

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

如何找到这些插值函数与水平线y = y0交叉的x值?以下是带有y0 = 2.85的图形化显示.

How can I find x-values where these interpolation functions cross a horizontal line y = y0? The following is a graphical illustration with 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)

我知道有关该主题的一些先前主题，例如

I am aware of a few previous threads on this topic, like

• 通过简单的拟合预测x值并在绘图中进行注释
• 使用拟合模型从Y值预测X值

• predict x values from simple fitting and annoting it in the plot
• Predict X value from Y value with a fitted model

建议我们简单地反转x和y，对(y, x)进行插值，然后计算在y = y0处的插值.

It is suggested that we simply reverse x and y, do an interpolation for (y, x) and compute the interpolated value at y = y0.

但是，这是一个虚假的想法.假设y = f(x)是(x, y)的插值函数，则此想法仅在f(x)是x的单调函数时有效，因此f是可逆的.否则x不是y的函数，并且对(y, x)进行插值是没有意义的.

However, this is a bogus idea. Let y = f(x) be an interpolation function for (x, y), this idea is only valid when f(x) is a monotonic function of x so that f is invertible. Otherwise x is not a function of y and interpolating (y, x) makes no sense.

用我的示例数据进行线性插值，这个假的想法给出了

Taking the linear interpolation with my example data, this fake idea gives

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

首先，根数不正确.我们从图中看到了两个根(在左侧)，但是代码仅返回一个.其次，它不是正确的根，

First of all, the number of roots is incorrect. We see two roots from the figure (on the left), but the code only returns one. Secondly, it is not a correct root, as

f1(fake_root)
#[1] 2.906103

不是2.85.

我已经在如何根据R中的roxfun()之后输入的y值估算x值，对此我做了首次尝试./a>.该解决方案对于线性插值来说是稳定的，但对于非线性插值来说不一定是稳定的.我现在正在寻找一种稳定的解决方案，特别是对于三次插值样条曲线.

I have made my first attempt on this general problem at How to estimate x value from y value input after approxfun() in R. The solution turns out stable for linear interpolation, but not necessarily stable for non-linear interpolation. I am now looking for a stable solution, specially for a cubic interpolation spline.

有时在单变量线性回归y ~ x或单变量非线性回归y ~ f(x)之后，我们要为目标y重新求解x.此问与答A是一个示例，并吸引了许多答案:解决最佳拟合多项式并绘制下拉线，但没有一个是真正的自适应性或易于在实践中使用.

Sometimes after a univariate linear regression y ~ x or a univariate non-linear regression y ~ f(x) we want to backsolve x for a target y. This Q & A is an example and has attracted many answers: Solve best fit polynomial and plot drop-down lines, but none is truly adaptive or easy to use in practice.

• 使用polyroot接受的答案仅适用于简单的多项式回归；
• 将二次公式用于解析解的答案仅适用于二次多项式；
• 我使用predict和uniroot所作的回答总体上可行，但并不方便，因为在实践中使用uniroot需要与用户进行交互(请参阅中的Uniroot解决方案，有关uniroot的更多信息.

• The accepted answer using polyroot only works for a simple polynomial regression;
• Answers using quadratic formula for an analytical solution only works for a quadratic polynomial;
• My answer using predict and uniroot works in general, but is not convenient, as in practice using uniroot needs interaction with users (see Uniroot solution in R for more on uniroot).

如果有一个自适应且易于使用的解决方案，那将是非常好的.

It would be really good if there is an adaptive and easy-to-use solution.

推荐答案

首先，让我复制.

## 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"方法返回的三次插值样条线，以下函数提供了稳定的数值解.

For cubic interpolation splines returned by stats::splinefun with methods "fmm", "natrual", "periodic" and "hyman", the following function provides a stable numerical solution.

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 的根.这是有效的，因为三次插值样条只是许多分段三次多项式.我在如何在R中保存和加载样条插值函数的答案?显示了如何获取分段多项式系数，因此使用polyroot简单明了.

It uses polyroot piecewise, first finding all roots on complex field, then retaining only real ones on the piecewise interval. This works because a cubic interpolation spline is just a number of piecewise cubic polynomials. My answer at How to save and load spline interpolation functions in R? has shown how to obtain piecewise polynomial coefficients, so using polyroot is straightforward.

使用问题中的示例数据，RootSpline1和RootSpline3都可以正确识别所有根.

Using the example data in the question, both RootSpline1 and RootSpline3 correctly identify all roots.

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

这篇关于在给定y值的情况下获得x值:线性/非线性插值函数的一般根查找的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！