2

I would like to use approxfun or a similar method to estimate a function for a curve and then estimate what the x value of a point would be given a known y value.

Here is a simplified example.

y <- seq(from=1, to =10, by = 1)

x <-seq(from=0.1, to =1, by = 0.1)

fun <- approxfun(x,y)

I can approximate a y value given a known x value with the following command:

fun(0.65)
#[1] 6.5

But how can I do the reverse, i.e., solve x from 6.5 = approxfun(x)?

Thanks for any help on my most likely dumb question.

Zheyuan Li
  • 71,365
  • 17
  • 180
  • 248
user005122
  • 23
  • 3

1 Answers1

1

Analytical solution for linear interpolation (stable)

Suppose we have some (x, y) data. After a linear interpolation find all x such that the value of the interpolant equals y0.

## with default value y0 = 0, it finds all roots of the interpolant
RootLinearInterpolant <- function (x, y, y0 = 0) {
  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)
  ## analytically root finding
  xk <- x[k] - z[k] * (x[k + 1] - x[k]) / (z[k + 1] - z[k])
  xk
  }

A more complicated example and test.

set.seed(0)
x <- sort(runif(10, 0, 10))
y <- rnorm(10, 3, 1)
y0 <- 2.5
xk <- RootLinearInterpolant(x, y, y0)
#[1] 3.375952 8.515571 9.057991

plot(x, y, "l"); abline(h = y0, lty = 2)
points(xk, rep.int(y0, length(xk)), pch = 19)

Numerical root finding for non-linear interpolation (not necessarily stable)

## suppose that f is an interpolation function of (x, y)
## this function finds all x, such that f(x) = y0
## with default value y0 = 0, it finds all roots of the interpolant
RootNonlinearInterpolant <- function (x, y, f, y0 = 0) {
  if (is.unsorted(x)) {
     ind <- order(x)
     x <- x[ind]; y <- y[ind]
     }
  z <- y - y0
  k <- which(z[-1] * z[-length(z)] < 0)
  nk <- length(k)
  xk <- numeric(nk)
  F <- function (x) f(x) - y0
  for (i in 1:nk) xk[i] <- uniroot(F, c(x[k[i]], x[k[i] + 1]))$root
  xk
  }

Try a natural cubic spline interpolation.

## cubic spline interpolation
f <- splinefun(x, y)
xk <- RootNonlinearInterpolant(x, y, f, y0)
#[1] 3.036643 8.953352 9.074306

curve(f, from = min(x), to = max(x))
abline(v = x, lty = 3)  ## signal pieces
abline(h = y0)
points(xk, rep.int(y0, length(xk)), pch = 20)

We see that that RootNonlinearInterpolant misses two crossover points on the 3rd piece.

RootNonlinearInterpolant relies on uniroot so the search is more restricted. Only if the sign of y - y0 changes on adjacent knots a uniroot is called. Clearly this does not hold on the 3rd piece. (Learn more about uniroot at Uniroot solution in R.)

Also note that uniroot only returns a single root. So the most stable situation is when the interpolant is monotone on the piece so a unique root exists. If there are actually multiple roots, uniroot would only find one of them.

Zheyuan Li
  • 71,365
  • 17
  • 180
  • 248