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:
1 2 3 4 5 6 7 8 9 10 11 12 |
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:
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:
1 2 3 4 5 6 |
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:
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].
1 2 3 4 5 6 7 8 9 10 11 |
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:
¿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 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.np.polyfit
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:
1 2 3 4 5 6 7 8 9 10 11 |
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:
1 2 |
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:
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:
1 2 3 |
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:
1 2 3 4 5 6 7 |
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í:
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
- RIVAS, Damián; VÁZQUEZ, Carlos. Elementos de Cálculo Numérico. ADI, 2010.
- BERRUT, Jean-Paul; TREFETHEN, Lloyd N. Barycentric lagrange interpolation. Siam Review, 2004, vol. 46, no 3, p. 501-517.
- 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]
muy bueno, esto me ha recordado la aproximación por mínimos cuadrados y cuando estudie la inversa de Moore-Penrose o g-inversa
Eso lo dejaremos para otro artículo 🙂 ¡Gracias por comentar!
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!
¡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 🙂
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 🙂
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! 🙂
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 🙂
En realidad, el secreto estaba en los colores básicamente: modifiqué el matplotlibrc para cambiar algunos de ellos y ahora los resultados me gustan mucho más 😛
https://gist.github.com/Juanlu001/6245967
¡Y ya no me quedan más ases en la manga! A este lado del Atlántico es hora ya de ir a dormir 😛 ¡Saludos!
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! 🙂
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!
¡Gracias a ti! Me alegro de que te haya resultado ameno, ¡un saludo! 😀
Reblogueó esto en Jack & Hack.
Pingback: Teoría de control en Python con SciPy (I) | Pybonacci
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
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!
Pingback: ¿Cómo borrar por encima de una línea en matplotlib? | Pybonacci
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.
¡Hola Eduardo, me alegro mucho! ¿Hablas de la ETSIN de la UPM? Si es así, somos vecinos 😉 ¡Un saludo!
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 🙁
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)
¿Donde consigo y como uso el modulo para hacer funcionar np.interp?
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
Que libreria tengo que importar para poder usar la funcion “np.interp”, porque no me funciona 🙁
Tienes que importar numpy, revisa el artículo. ¡Un saludo!
Muchas gracias por vuestros comentarios, para más dudas os recomendamos http://es.stackoverflow.com/. ¡Un saludo! 🙂
Los comentarios están cerrados.