El ajuste de curvas es el proceso de construir una curva (función), que sea el mejor ajuste a una serie de puntos. Las curvas ajustadas pueden ser usadas como asistencia en la visualización de datos, para inferir valores de una función donde no hay datos disponibles, y para resumir la relación entre variables.
Referencia:
Consideremos un polinomio de grado uno:
$$y = \beta_1 x + \beta_0.$$Esta es una línea recta que tiene pendiente a_1. Sabemos que habrá una línea conectando dos puntos cualesquiera. Por tanto, una ecuación polinómica de primer grado es un ajuste perfecto entre dos puntos.
Si consideramos ahora un polinomio de segundo grado,
$$y = \beta_2 x^2 + \beta_1 x + \beta_0,$$este se ajustará exactamente a tres puntos. Si aumentamos el grado de la función a la de un polinomio de tercer grado, obtenemos:
$$y = \beta_3 x^3 + \beta_2 x^2 + \beta_1 x + \beta_0,$$que se ajustará a cuatro puntos.
Ejemplos
Solución
In [1]:
import numpy as np
import matplotlib.pyplot as plt
In [2]:
P1 = [0, 1]
P2 = [1, 0]
X = np.array([[1, 0], [1, 1]])
y = np.array([1, 0])
b0, b1 = np.linalg.inv(X).dot(y)
b0, b1
Out[2]:
In [3]:
x = np.linspace(-0.2, 1.2, 100)
y = b1*x+b0
plt.figure(figsize=(8,6))
plt.plot([P1[0], P2[0]], [P1[1], P2[1]], 'r*', label = 'puntos')
plt.plot(x, y, 'b', label = 'recta ajustada')
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.legend(loc = 'best')
plt.grid()
plt.show()
In [4]:
P1 = [-1, 1]
P2 = [0, 0]
P3 = [1, 1]
X = np.array([[1, -1, 1], [1, 0, 0], [1, 1, 1]])
y = np.array([1, 0, 1])
b0, b1, b2 = np.linalg.inv(X).dot(y)
b0, b1, b2
Out[4]:
In [5]:
x = np.linspace(-1.2, 1.2, 100)
y = b2*x**2+b1*x+b0
plt.figure(figsize=(8,6))
plt.plot([P1[0], P2[0], P3[0]], [P1[1], P2[1], P3[1]], 'r*', label = 'puntos')
plt.plot(x, y, 'b', label = 'parábola ajustada')
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.legend(loc = 'best')
plt.grid()
plt.show()
Las curvas están completamente determinadas por los puntos (datos limpios, suficientes y necesarios).
Esto se traduce en que, al llevar el problema a un sistema de ecuaciones lineales, existe una única solución: no hay necesidad, ni se puede optimizar nada.
¿Tendremos datos así de 'bonitos' en la vida real?
La realidad es que los datos que encontraremos en nuestra vida profesional se parecen más a esto...
In [37]:
x = np.linspace(0, 1, 30)
y = 10*x + 2 + np.random.randn(30)
plt.figure(figsize=(8,6))
plt.plot(x, y, '*b')
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.grid()
plt.show()
Consideramos que tenemos un conjunto de n pares ordenados de datos $(x_i,y_i)$, para $i=1,2,3,\dots,n$.
Consideramos entonces ajustes de la forma $\hat{f}(x) = \beta_0+\beta_1 x = \left[1 \quad x\right]\left[\begin{array}{c} \beta_0 \\ \beta_1 \end{array}\right]=\left[1 \quad x\right]\boldsymbol{\beta}$ (lineas rectas).
Para decir 'mejor', tenemos que definir algún sentido en que una recta se ajuste mejor que otra.
Mínimos cuadrados: el objetivo es seleccionar los coeficientes $\boldsymbol{\beta}=\left[\beta_0 \quad \beta_1 \right]^T$, de forma que la función evaluada en los puntos $x_i$ ($\hat{f}(x_i)$) aproxime los valores correspondientes $y_i$.
La formulación por mínimos cuadrados, encuentra los $\boldsymbol{\beta}=\left[\beta_0 \quad \beta_1 \right]^T$ que minimiza $$\sum_{i=1}^{n}(y_i-\hat{f}(x_i))^2=\sum_{i=1}^{n}(y_i-\left[1 \quad x_i\right]\boldsymbol{\beta})^2=\left|\left|\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right|\right|^2,$$
donde $\boldsymbol{y}=\left[y_1\quad\dots\quad y_n\right]$, y $\boldsymbol{X}=\left[\begin{array}{ccc}1 & x_1\\ \vdots & \vdots \\ 1 & x_n\end{array}\right].$ Esto es,
$$\boldsymbol{\beta}^{ls} = \arg \min_{\boldsymbol{\beta}} \left|\left|\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right|\right|^2$$Para llevar a cabo la anterior minimización, el módulo pyomo_utilities.py ha sido actualizado con la función curve_polyfit(x, y, order, reg = None, robust = False). Esta función recibe:
Para los datos anteriores queríamos ajustar una recta (order = 1)...
In [2]:
import pyomo_utilities
In [39]:
beta = pyomo_utilities.curve_polyfit(x, y, 1)
In [40]:
yhat = beta[0]+beta[1]*x
plt.figure(figsize=(8,6))
plt.plot(x, y, '*b', label = 'datos')
plt.plot(x, yhat, '-r', label = 'ajuste')
plt.legend(loc = 'best')
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.grid()
plt.show()
In [41]:
beta
Out[41]:
Note que la pendiente es aproximadamente $10$ y el intercepto es aproximadamente $2$.
La anterior idea se puede extender a ajuste polinomial...
In [42]:
n = 100
x = np.linspace(np.pi/3, 5*np.pi/3, n)
y = 4*np.sin(x) + 0.5*np.random.randn(n)
plt.figure(figsize=(8,6))
plt.plot(x, y, '*b')
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.grid()
plt.show()
¿Se ajustará bien una recta?
In [43]:
beta1 = pyomo_utilities.curve_polyfit(x, y, 1)
yhat1 = beta1[0]+beta1[1]*x
plt.figure(figsize=(8,6))
plt.plot(x, y, '*b', label = 'datos')
plt.plot(x, yhat1, '-r', label = 'ajuste 1')
plt.legend(loc = 'best')
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.grid()
plt.show()
In [44]:
beta1
Out[44]:
Creo que no. Pero, ¿Y una parábola?
In [45]:
beta2 = pyomo_utilities.curve_polyfit(x, y, 2)
yhat2 = beta2[0]+beta2[1]*x+beta2[2]*x**2
plt.figure(figsize=(8,6))
plt.plot(x, y, '*b', label = 'datos')
plt.plot(x, yhat1, '-r', label = 'ajuste 1')
plt.plot(x, yhat2, '-g', label = 'ajuste 2')
plt.legend(loc = 'best')
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.grid()
plt.show()
In [46]:
beta2
Out[46]:
Tampoco. Quizá un polinomio cúbico...
In [47]:
beta3 = pyomo_utilities.curve_polyfit(x, y, 3)
yhat3 = beta3[0]+beta3[1]*x+beta3[2]*x**2+beta3[3]*x**3
plt.figure(figsize=(8,6))
plt.plot(x, y, '*b', label = 'datos')
plt.plot(x, yhat1, '-r', label = 'ajuste 1')
plt.plot(x, yhat2, '-g', label = 'ajuste 2')
plt.plot(x, yhat3, '-k', label = 'ajuste 3')
plt.legend(loc = 'best')
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.grid()
plt.show()
In [48]:
beta3
Out[48]:
Mucho mejor. Entonces, ¿mientras más se suba el orden mejor la aproximación?
¡Cuidado! OVERFITTING...
In [49]:
beta20 = pyomo_utilities.curve_polyfit(x, y, 20)
yhat20 = np.array([x**i for i in range(21)]).T.dot(beta20)
plt.figure(figsize=(8,6))
plt.plot(x, y, '*b', label = 'datos')
plt.plot(x, yhat1, '-r', label = 'ajuste 1')
plt.plot(x, yhat2, '-g', label = 'ajuste 2')
plt.plot(x, yhat3, '-k', label = 'ajuste 3')
plt.plot(x, yhat20, '-c', label = 'ajuste 20')
plt.legend(loc = 'best')
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.grid()
plt.show()
In [50]:
beta20
Out[50]:
¡Cuidado! ver el tamaño de algunos coeficientes. Cuando los coeficientes son grandes, ¿qué pasa?
Es conveniente ver el error como función del orden del polinomio... selección de modelos
In [51]:
e_ms = []
for i in range(10):
beta = pyomo_utilities.curve_polyfit(x, y, i+1);
yhat = np.array([x**j for j in range(i+2)]).T.dot(beta)
e_ms.append(sum((y - yhat)**2))
plt.figure(figsize=(8,5))
plt.plot(np.arange(10)+1, e_ms, 'o')
plt.xlabel('orden')
plt.ylabel('error')
plt.show()
En efecto, parece que con $3$ es suficiente.
Vimos que la solución de mínimos cuadrados es: $$\boldsymbol{\beta}^{ls} = \arg \min_{\boldsymbol{\beta}} \left|\left|\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right|\right|^2.$$
Sin embargo, si crecemos el orden del modelo hay overfitting y algunos coeficientes óptimos $\boldsymbol{\beta}$ crecen muchísimo. Que un coeficiente sea muy grande, significa que se le da mucha importancia a alguna característica (que quizá sea ruido... no sirve para predecir).
La regularización consiste en penalizar la magnitud de los coeficientes $\boldsymbol{\beta}$ en el problema de optimización, para que no crezcan tanto.
In [52]:
beta20_ridge = pyomo_utilities.curve_polyfit(x, y, 20, reg_mode = 'ridge', reg_coef = 0.1)
yhat20_ridge = np.array([x**i for i in range(21)]).T.dot(beta20_ridge)
plt.figure(figsize=(8,6))
plt.plot(x, y, '*b', label = 'datos')
plt.plot(x, yhat20, '-c', label = 'ajuste 20')
plt.plot(x, yhat20_ridge, '--r', label = 'ajuste 20_ridge')
plt.legend(loc = 'best')
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.grid()
plt.show()
In [53]:
beta20_ridge
Out[53]:
La norma 1 no es más que la suma de los valores absolutos de las componentes $\left|\left|\boldsymbol{\beta}\right|\right|_1=\sum_{j=0}^m\left|\beta_j\right|$.
In [54]:
beta20_lasso = pyomo_utilities.curve_polyfit(x, y, 20, reg_mode = 'lasso', reg_coef = 0.001)
yhat20_lasso = np.array([x**i for i in range(21)]).T.dot(beta20_lasso)
plt.figure(figsize=(8,6))
plt.plot(x, y, '*b', label = 'datos')
plt.plot(x, yhat20, '-c', label = 'ajuste 20')
plt.plot(x, yhat20_ridge, '--r', label = 'ajuste 20_ridge')
plt.plot(x, yhat20_lasso, '--k', label = 'ajuste 20_lasso')
plt.legend(loc = 'best')
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.grid()
plt.show()
In [55]:
beta20_lasso
Out[55]:
In [3]:
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(0, 1, 30)
y = 10*x + 2 + np.random.randn(30)
y[0] = 16
y[-1] = 0
plt.figure(figsize=(8,6))
plt.plot(x, y, '*b')
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.grid()
plt.show()
In [4]:
beta = pyomo_utilities.curve_polyfit(x, y, 1)
yhat = beta[0]+beta[1]*x
plt.figure(figsize=(8,6))
plt.plot(x, y, '*b', label = 'datos')
plt.plot(x, yhat, '-r', label = 'ajuste')
plt.legend(loc = 'best')
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.grid()
plt.show()
In [5]:
beta
Out[5]:
Si estos puntos que parecen ser atípicos, hacen parte de una 'mala medición', vemos que el ajuste que obtenemos a los otros puntos es muy pobre...
¿Cómo podemos evitar esto? La respuesta es ajuste robusto.
In [6]:
beta = pyomo_utilities.curve_polyfit(x, y, 1, robust=True)
yhat = beta[0]+beta[1]*x
plt.figure(figsize=(8,6))
plt.plot(x, y, '*b', label = 'datos')
plt.plot(x, yhat, '-r', label = 'ajuste')
plt.legend(loc = 'best')
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.grid()
plt.show()
In [7]:
beta
Out[7]:
Mejor...
Abrir un nuevo notebook, llamado Tarea7_ApellidoNombre y subirlo a moodle en el espacio habilitado. Si no terminan esto en clase, tienen hasta mañana a las 23:00.
In [30]:
def f(x):
return np.exp(-x**2/2)/np.sqrt(2*np.pi)
In [31]:
x = np.linspace(-3, 3)
y = f(x)
plt.figure(figsize=(8,6))
plt.plot(x, y, '*b', label = 'datos')
plt.legend(loc = 'best')
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.grid()
plt.show()
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]: