obtener el valor de x dado el valor de y: búsqueda raíz general para la función de interpolación lineal/no lineal

regression interpolation

Question

Estoy interesado en un problema general de búsqueda de raíz para una función de interpolación.

Supongamos que tengo los siguientes datos (x, y) :

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

así como una interpolación lineal y una interpolación spline cúbica:

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

¿Cómo puedo encontrar los valores de x donde estas funciones de interpolación cruzan una línea horizontal y = y0 ? La siguiente es una ilustración gráfica con 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)

Soy consciente de algunos hilos anteriores sobre este tema, como

• Predecir los valores de x a partir de un ajuste simple y anotarlos en la gráfica.
• Predecir el valor X a partir del valor Y con un modelo ajustado

Se sugiere que simplemente invirtamos x e y , hagamos una interpolación para (y, x) y calculamos el valor interpolado en y = y0 .

Sin embargo, esta es una idea falsa. Sea y = f(x) una función de interpolación para (x, y) , esta idea solo es válida cuando f(x) es una función monotónica de x para que f sea ​​invertible. De lo contrario, x no es una función de y y la interpolación (y, x) no tiene sentido.

Tomando la interpolación lineal con mis datos de ejemplo, esta idea falsa da

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

En primer lugar, el número de raíces es incorrecto. Vemos dos raíces en la figura (a la izquierda), pero el código solo devuelve una. En segundo lugar, no es una raíz correcta, como

f1(fake_root)
#[1] 2.906103

No es 2.85.

Hice mi primer intento en este problema general en Cómo estimar el valor de x a partir de la entrada del valor y después de aproxfun () en R. La solución resulta estable para la interpolación lineal, pero no necesariamente estable para la interpolación no lineal. Ahora estoy buscando una solución estable, especialmente para una spline de interpolación cúbica.

¿Cómo puede una solución ser útil en la práctica?

A veces, después de una regresión lineal univariada y ~ x xo una regresión univariable no lineal y ~ f(x) , queremos dar una vuelta atrás a x para un objetivo y . Este Q & A es un ejemplo y ha atraído muchas respuestas: resuelva el polinomio de mejor ajuste y dibuje líneas desplegables , pero ninguna es realmente adaptable o fácil de usar en la práctica.

• La respuesta aceptada utilizando polyroot solo funciona para una regresión polinomial simple;
• Las respuestas que usan la fórmula cuadrática para una solución analítica solo funcionan para un polinomio cuadrático;
• Mi respuesta con predict y uniroot funciona en general, pero no es conveniente, ya que en la práctica uniroot necesario interactuar con los usuarios (ver la solución Uniroot en R para más información sobre uniroot ).

Sería realmente bueno si hubiera una solución adaptable y fácil de usar.

Answer 1

Dados los puntos de datos y la función de spline como anteriormente, simplemente aplique findzeros() desde el paquete pracma .

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

Answer 2

En primer lugar, permítame copiar en la solución estable para la interpolación lineal propuesta en mi respuesta anterior .

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

Para las splines de interpolación cúbicas devueltas por stats::splinefun con los métodos "fmm" , "natrual" , "periodic" y "hyman" , la siguiente función proporciona una solución numérica estable.

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
}

Utiliza la polyroot partes, primero encuentra todas las raíces en un campo complejo y luego retiene solo las reales en el intervalo por partes. Esto funciona porque una spline de interpolación cúbica es solo un número de polinomios cúbicos por partes. Mi respuesta en ¿Cómo guardar y cargar funciones de interpolación spline en R? ha demostrado cómo obtener coeficientes polinomiales por partes, por lo que usar polyroot es sencillo.

Usando los datos de ejemplo en la pregunta, tanto RootSpline1 como RootSpline3 identifican correctamente todas las raíces.

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