Search code examples
rregressioninterpolationspline

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

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

Suppose I have the following (x, y) data:

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

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

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

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

is not 2.85.

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.

How can a solution be useful in practice?

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.

  • 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.

Solution

  • First of all, let me copy in the stable solution for linear interpolation proposed in my previous answer.

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

    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
  }

    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.

    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