Ajuste e interpolación unidimensionales básicos en Python con SciPy

Introducción

En este artículo vamos a ver una introducción a cómo hacer ajustes e interpolaciones en Python utilizando NumPy y los módulos interpolate y optimize de SciPy.

Ajustes de curvas e interpolaciones son dos tareas básicas que realizaremos con mucha frecuencia. Por ejemplo, cuando recojamos los datos de un experimento: sabemos que se tienen que comportar como una parábola, pero obviamente por errores de medición u otro tipo no obtenemos una parábola exactamente. En este caso necesitaremos realizar un ajuste de los datos, conocido el modelo (una curva de segundo grado en este caso).

En otras ocasiones dispondremos de una serie de puntos y querremos construir una curva que pase por todos ellos. En este caso lo que queremos es realizar una interpolación: si tenemos pocos puntos podremos usar un polinomio, y en caso contrario habrá que usar trazadores (splines en inglés). Vamos a empezar por este último método.

Si deseas consultar el código completo (incluyendo el que genera las figuras) puedes ver el notebook que usé para redactar el artículo.

En esta entrada se han usado python 3.3.2, numpy 1.7.1, scipy 0.12.0 y matplotlib 1.3.0.

Interpolación

Polinomios no, ¡gracias!

Lo primero que vamos a hacer va a ser desterrar la idea de que, sea cual sea el número de puntos que tengamos, podemos construir un polinomio que pase por todos ellos «y que lo haga bien». Si tenemos \(N\) puntos nuestro polinomio tendrá que ser de grado menor o igual que \(N - 1\), pero cuando \(N\) empieza a ser grande (del orden de 10 o más) a menos que los puntos estén muy cuidadosamente elegidos el polinomio oscilará salvajemente. Esto se conoce como fenómeno de Runge.

Para ver esto podemos estudiar el clásico ejemplo que dio Runge: tenemos la función

\(\displaystyle f(x) = \frac{1}{1 + x^2}\)

veamos qué sucede si la interpolamos en nodos equiespaciados. Para ello vamos a usar la función barycentric_interpolate (según Berrut y Trefethen [II] «[El método de interpolación baricéntrica] merece ser conocido como el método estándar de interpolación polinómica»). Esta función recibe tres argumentos:

  • una lista de coordenadas x_i de los nodos,
  • una lista de coordenadas y_i de los nodos, y
  • un array x donde evaluar el polinomio interpolante que resulta.

El código será este:

import numpy as np
from scipy.interpolate import barycentric_interpolate
def runge(x):
    """Función de Runge."""
    return 1 / (1 + x ** 2)
N = 11  # Nodos de interpolación
xp = np.arange(11) - 5  # -5, -4, -3, ..., 3, 4, 5
fp = runge(xp)
x = np.linspace(-5, 5)
y = barycentric_interpolate(xp, fp, x)

Y este es el resultado:

Fenómeno de Runge
Con 11 nodos equiespaciados en este ejemplo dado por Runge el polinomio interpolante diverge en los extremos.

Y no os quiero contar nada si escogemos 20 o 100 puntos.

Existe una forma de mitigar este problema, que es, como ya hemos dicho, «escogiendo los puntos cuidadosamente». Una de las formas es elegir las raíces de los polinomios de Chebyshev, que podemos construir en NumPy usando el módulo polynomial.chebyshev. Por ejemplo, si queremos como antes 11 nodos tendremos que escoger el polinomio de Chebyshev de grado 11:

from numpy.polynomial import chebyshev
coeffs_cheb = [0] * 11 + [1]  # Solo queremos el elemento 11 de la serie
T11 = chebyshev.Chebyshev(coeffs_cheb, [-5, 5])
xp_ch = T11.roots()
# -4.949, -4.548, -3.779, -2.703, ..., 4.548, 4.949

Utilizando estos puntos, la cosa no queda tan mal:

Comparación
Usando nodos de Chebyshev la interpolación es mucho mejor.

Aun así, aún tenemos varios problemas:

  • El polinomio sigue oscilando, y esto puede no ser deseable.
  • No siempre podemos escoger los puntos como nosotros queramos.

Por tanto, desde ya vamos a abandonar la idea de usar polinomios y vamos a hablar de trazadores (splines en inglés).

Trazadores

Los trazadores o splines no son más que curvas polinómicas definidas a trozos, normalmente de grado 3 (casi nunca mayor de 5). Al ser cada uno de los trozos de grado pequeño se evita el fenómeno de Runge, y si se «empalman» los trozos inteligentemente la curva resultante será suave (matemáticamente: diferenciable) hasta cierto punto. Cuando queremos una curva que pase por todos los puntos disponibles un trazador es justamente lo que necesitamos.

El trazador más elemental, el lineal (grado 1), se puede construir rápidamente en NumPy usando np.interp. El más común, el trazador cúbico (grado 3) se puede construir con la clase scipy.interpolate.InterpolatedUnivariateSpline. Si pasamos a esta clase un argumento k podemos especificar el grado del trazador (entre 1 y 5). Como ejemplo vamos a tomar los datos de la silueta del pato de Villafuerte [III].

from scipy.interpolate import InterpolatedUnivariateSpline
# Pato
P = [(0.9, 1.3), (1.3, 1.5), (1.9, 1.8), (2.1,2.1), (2.6, 2.6), (3.0, 2.7),
     (3.9, 2.3), (4.4, 2.1), (4.8, 2.0), (5.0, 2.1), (6, 2.2), (7, 2.3),
     (8, 2.2), (9.1, 1.9), (10.5, 1.4), (11.2, 0.9), (11.6, 0.8), (12, 0.6),
     (12.6, 0.5), (13, 0.4), (13.2, 0.2)]
xi, yi = zip(*P)  # 21 puntos de interpolación
x = np.linspace(min(xi), max(xi), num=1001)  # Dominio
y1d = np.interp(x, xi, yi)
#y1d = InterpolatedUnivariateSpline(xi, yi, k=1)(x)  # Mismo resultado
ysp = InterpolatedUnivariateSpline(xi, yi)(x)  # Llamamos a la clase con x

Nota: ¿Quieres saber el truco de zip(*P)? 😉

Y si representamos el resultado obtenemos esto:

Trazadores lineal y cúbico.
Trazadores lineal y cúbico.

¿Alguien se anima a enviarnos una gráfica de cómo quedaría la interpolación si usásemos un polinomio de grado 20? 😉

En ocasiones, sin embargo, puede que no necesitemos un trazador que pase por todos los puntos, sino una curva o un modelo más sencillo que aproxime una serie de puntos, tratando de cometer el mínimo error posible. Si quieres saber cómo hacer esto, ¡sigue leyendo!

Ajuste de curvas

Ajuste polinómico

El ajuste más básico en el que podemos pensar es el ajuste polinómico: buscamos un polinomio que aproxime los datos con el menor error posible. Para ello utilizaremos la función polynomial.polyfit del paquete polynomial de NumPy.

Nota: La función np.polyfit es diferente a la que vamos a usar aquí y está obsoleta, aparte de que tiene el convenio contrario para los coeficientes. Se recomienda no usarla. Ya sé que la otra tiene un nombre un poco largo y que los ejemplos de la documentación tienen fallos, pero es lo que hay.

Esta función recibe tres argumentos obligatorios:

  • una lista de coordenadas x de los puntos,
  • una lista de coordenadas y de los puntos, y
  • el grado deg del polinomio interpolante.

Vamos a ver un ejemplo real con el que me encontré hace unos meses: hallar la polar parabólica aproximada de un avión. Para ello podéis utilizar estos datos.

La polar de un avión es la relación entre la sustentación y la resistencia aerodinámica del mismo. Su forma teórica es:

\(\displaystyle C_D = C_{D0} + k C_L^2\)

podríamos estar tentados de intentar un ajuste parabólico, pero vemos que en realidad no aparece término lineal. Si llamamos \(y = C_D\) y \(x = C_L^2\) tenemos:

\(\displaystyle y = y_0 + k x\)

con lo que podemos realizar un ajuste lineal. Por otro lado, tengo que descartar los puntos que están más allá de la condición de entrada en pérdida (después del máximo del coeficiente de sustentación), porque esos no cuadran con el modelo teórico. Este es el código:

import numpy.polynomial as P
# Cargamos los datos
data = np.loadtxt("polar.dat")
_, C_L, C_D = data
# Descarto los datos que no me sirven
stall_idx = np.argmax(C_L)
y = C_D[:stall_idx + 1]
x = C_L[:stall_idx + 1] ** 2
# Ajuste lineal, devuelve los coeficientes en orden creciente
C_D0, k = P.polynomial.polyfit(x, y, deg=1)
print(C_D0, k)

Una vez hemos obtenido los dos coeficientes, no hay más que evaluar el polinomio en un cierto dominio usando la función polynomial.polyval, que acepta como argumentos

  • el dominio donde queremos evaluar la función y
  • una lista de coeficientes de grado creciente, tal y como los devuelve polyfit.

El código es simplemente:

C_L_dom = np.linspace(C_L[0], C_L[stall_idx], num=1001)
C_D_int = P.polynomial.polyval(C_L_dom ** 2, (C_D0, k))

Y este es el resultado que obtenemos:

Polar parabólica y datos reales.
Polar parabólica y datos reales.

En la figura se aprecia perfectamente cómo he descartado los puntos más allá del máximo y cómo la parábola, aun no pasando por todos los puntos (tal vez no pase por ninguno) aproxima bastante bien los datos que tenemos. ¡Bien!

General

En ocasiones las cosas son más complicadas que un polinomio (sí lectores, así es la vida). Pero no pasa nada, porque con la función scipy.optimize.curve_fit podemos ajustar una serie de datos a cualquier modelo que se nos ocurra, no importa qué tan complicado sea. Sin ir más lejos, tomando el ejemplo de la documentación, supongamos que tenemos unos datos que se ajustan al modelo

\(\displaystyle A e^{-B x^2} + C\)

en Python nuestro modelo será una función que recibirá como primer argumento x y el resto serán los parámetros del mismo:

def func(x, A, B, C):
    """Modelo para nuestros datos."""
    return A * np.exp(-B * x ** 2) + C

Ahora solo necesitamos algunos datos (añadiremos un poco de ruido gaussiano para que tenga más gracia) y podemos probar el ajuste. La función scipy.optimize.curve_fit recibe como argumentos:

  • el modelo func para los datos,
  • una lista de coordenadas xdata de los puntos, y
  • una lista de coordenadas ydata de los puntos.

Así realizamos el ajuste:

from scipy.optimize import curve_fit
Adat, Bdat, Cdat = 2.5, 1.3, 0.5
xdat = np.linspace(-2, 4, 12)
ydat = func(xdat, Adat, Bdat, Cdat) + 0.2 * np.random.normal(size=len(xdat))
(A, B, C), _ = curve_fit(func, xdat, ydat)
print(A, B, C)

Y el resultado queda así:

Hemos especificado nuestro modelo, y este es el ajuste.
Hemos especificado nuestro modelo, y este es el ajuste.

Fácil, ¿no?

Mínimos cuadrados

Todas estas funciones emplean la solución por mínimos cuadrados de un sistema lineal. Nosotros podemos acceder a esta solución utilizando la función scipy.optimize.leastsq, pero como es más general y este artículo ya se ha extendido bastante lo vamos a dejar aquí, de momento 😉

Y tú, ¿te animas ya a realizar ajustes e interpolaciones con Python? ¿Qué dificultades ves? ¿Cómo piensas que podríamos mejorar el artículo? ¡Cuéntanoslo en los comentarios! 🙂

Referencias

  1. RIVAS, Damián; VÁZQUEZ, Carlos. Elementos de Cálculo Numérico. ADI, 2010.
  2. BERRUT, Jean-Paul; TREFETHEN, Lloyd N. Barycentric lagrange interpolation. Siam Review, 2004, vol. 46, no 3, p. 501-517.
  3. VILLAFUERTE, Héctor F. Guías para Métodos Numéricos, parte 2 [en línea]. 2010. Disponible en web: <http://uvgmm2010.wordpress.com/guias/>. [Consulta: 15 de agosto de 2013]

Juan Luis Cano

Estudiante de ingeniería aeronáutica y con pasión por la programación y el software libre. Obsesionado con mejorar los pequeños detalles y con ganas de cambiar el mundo. Divulgando Python en español a través de Pybonacci y la asociación Python España.

More Posts - Website

Follow Me:
TwitterLinkedIn

25 thoughts on “Ajuste e interpolación unidimensionales básicos en Python con SciPy

  1. muy bueno, esto me ha recordado la aproximación por mínimos cuadrados y cuando estudie la inversa de Moore-Penrose o g-inversa

  2. Bueno, lo estoy leyendo y todavía no terminé, pero me ganó la ansiedad en felicitar y decir que está todo muy interesante!

    No entendí lo del zip(*P)… lo probé y funciona, pero no entiendo el asterisco… también me hubiese gustado ver los comandos para hacer los gráficos (con sus tipos de curvas, anotaciones, etc…), o por lo menos para el primero 🙂

    Abrazos, muchas gracias!

    1. ¡Hola Germán! Muchas gracias por comentar 😀 Tienes razón en las dos cosas: en que falta el código de los gráficos y que lo de zip(*P) es magia negra 😛

      He enlazado ya los gráficos aquí: http://nbviewer.ipython.org/6245476

      En cuanto a lo segundo, cuando a una función le pasas un argumento iterable con un asterisco delante, se «desempaqueta». Por ejemplo:

      [sourcecode language="python"]
      def foo(a, b):
      return a + b

      #foo([1, 2]) # Error! Falta un argumento
      foo(*[1, 2]) # 3

      [/sourcecode]

      Lo que haces con zip(*P), es equivalente a hacer zip((0.9, 1.3), (1.3, 1.5), …). Si no pones asterisco, queda zip([(0.9, 1.3), (1.3, 1.5), …]), que parece lo mismo pero ya ves que no es igual 😉

      Quería darte una explicación rápida. En inglés se llama “argument unpacking”. Seguro que si preguntas en Python Majibu (http://python.majibu.org/) alguno de los maestros que contestan ahí te darán una descripción excelente 🙂

      1. Genial ese notebook Juan… ya mismo me lo estoy bajando!

        Tu explicación del asterisco está muy clara, por más rápida que sea jeje… ya mismo me la apunto!

        Abrazos y gracias nuevamente 🙂

      2. Las gráficas casi seguro que no te van a quedar igual, porque aún tengo otro truco reservado… 😛 ¡Si te interesa saberlo coméntanoslo por aquí!

        ¡Un saludo y gracias a ti! 🙂

      3. Los dos primeros gráficos me quedaron casi iguales (todavía no probé otros), a no ser por las vocales acentuadas y los límites en el eje ‘y’… fue copy/paste del notebook hacia ipython… ahora me dejaste curioso jaja!

        Espero ver esos ases que tenés guardados en la manga 🙂

  3. Justamente te iba a decir en la respuesta anterior que soy daltónico por si tenía que ver con los colores, pero me olvidé… jajajaa… Gracias por el gist, lo voy a ver… buenas noches! 🙂

  4. Muchas gracias por el artículo, me ha gustado leerlo, y lo tendré en cuenta, creo que lo necesitaré en cosa de días 🙂 Muchas gracias por el notebook!

  5. Tengo un articulo que resulta dificil de entender pues esta en ingles y no soy muy bueno con los idiomas, es sobre interpolacion y se explica el metodo de lagrange, newton y muchos mas pero no se ha que se refiere, no se si pudieras ayudarme

    1. Hola, lo siento pero no puedo ayudarte con eso 🙁 De todos modos en las referencias que damos al final hay textos muy buenos en español y fácilmente accesibles. ¡Un saludo!

  6. Magnífico. En mi escuela (navales) se ha impartido un curso de Python para acercar al alumnado a este software, y tu artículo ha sido de gran ayuda para el entregable práctico que debíamos hacer. Muchas gracias.

    1. ¡Hola Eduardo, me alegro mucho! ¿Hablas de la ETSIN de la UPM? Si es así, somos vecinos 😉 ¡Un saludo!

  7. Hola, buena informacion. Una consulta cuando ejecuto el programa de ajuste polinomico, me sale que np no esta definido.
    Pd. Algunos de los metodos esta relacionado con el goodness of fit? Tengo un txt con datos en x e y; y quiero ajustar la recta, porque no sale una recta perfecta. Tengo 4 opciones, lineal, cuadractica, cubica y exponencial. Alguna idea? Soy nuevo en python 🙁

    1. Has escrito esto antes?
      import numpy as np

      Python funciona con modulos, si queremos usar las funciones de alguno de ellos primero tenemos que importarlo.

      (no estoy en un teclado espanol)

  8. Una pregunta alguien sabe como crear un funcion por teclado … algo asi como lo hace el matlab
    ejemplo :
    f(x)=x^2 +1
    evaluar f(3)
    pero debes ingresar por teclado “x**2+1”.
    mas especifico cuando uso el matlab
    puedo usar
    g=input(“ingrese funcion:”)
    f=inline(g,’x’)
    esto hara que f dependa x eso quiero hacer de la misma manera en python alguna idea … gracias

  9. Que libreria tengo que importar para poder usar la funcion “np.interp”, porque no me funciona 🙁

Comments are closed.