El tema principal en este módulo fue optimización. Al finalizar este módulo, se espera que ustedes tengan las siguientes competencias
- Realizar optimizaciones de funciones escalares en un dominio dado usando sympy.
- Dado un problema de programación lineal, llevarlo a la forma que vimos en clase y resolverlo.
- Ajustar curvas a conjuntos de puntos dados.
- Diseñar clasificadores binarios con regresión logística para conjuntos de datos linealmente separables.
En clase vimos cómo optimizar funciones escalares dado un invervalo cerrado finito utilizando sympy. Ustedes, además, realizaron una tarea donde hicieron una función genérica para optimizar cualquier función dada.
Recordamos en este ejemplo como optimizar este tipo de funciones.
Obtener el máximo y el mínimo absoluto de la función
$$f(x)=2x^4-16x^3+32x^2+5$$en el intervalo $[-1, 4.5]$. Graficar la función en este intervalo, junto con el punto donde ocurre el máximo (en color rojo) y el punto donde ocurre el mínimo (en color azul).
Les recuerdo nada más como imprimir en formato LaTeX.
In [1]:
# Librería de cálculo simbólico
import sympy as sym
# Para imprimir en formato TeX
from sympy import init_printing; init_printing(use_latex='mathjax')
¿Qué más librerías necesitamos?
In [2]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
¿Qué más sigue?
In [3]:
def f(x):
return 2*x**4-16*x**3+32*x**2+5
In [4]:
sym.var('x', real = True)
Out[4]:
In [5]:
f(x)
Out[5]:
In [6]:
df = sym.diff(f(x), x)
xc = sym.solve(df, x)
xc
Out[6]:
In [7]:
f(-1), f(4.5), f(xc[0]), f(xc[1]), f(xc[2])
Out[7]:
In [8]:
xnum = np.linspace(-1, 4.5, 100)
plt.figure(figsize=(8,6))
plt.plot(xnum, f(xnum), 'k', lw = 2, label = '$y=f(x)$')
plt.plot([-1], [f(-1)], '*r', ms = 10, label = '$max_{-1\leq x\leq 4.5}f(x)=55$')
plt.plot([xc[0], xc[2]], [f(xc[0]), f(xc[2])], '*b', ms = 10, label = '$min_{-1\leq x\leq 4.5}f(x)=5$')
plt.legend(loc='best')
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.show()
In [9]:
def A(theta):
return (1+sym.sin(theta))*sym.cos(theta)
In [10]:
sym.var('theta', real = True)
Out[10]:
In [11]:
A(theta)
Out[11]:
In [12]:
dA = sym.diff(A(theta), theta)
dA
Out[12]:
In [13]:
thetac = sym.solve(dA, theta)
thetac
Out[13]:
In [14]:
A(0), A(thetac[1]), A(np.pi/2)
Out[14]:
In [15]:
thetanum = np.linspace(0, np.pi/2, 100)
Anum = sym.lambdify([theta], A(theta), 'numpy')
plt.figure(figsize=(8,6))
plt.plot(thetanum, Anum(thetanum), 'k', lw = 2, label = '$y=A(\theta)$')
plt.plot([thetac[1]], [A(thetac[1])], '*r', ms = 10, label = '$max_{0\leq x\leq \pi/2}A(\theta)$')
plt.legend(loc='best')
plt.xlabel('$\theta$')
plt.ylabel('$y$')
plt.show()
En clase vimos cómo llevar los problemas de programación lineal a la forma \begin{equation} \begin{array}{ll} \min_{\boldsymbol{x}} & \boldsymbol{f}^T\boldsymbol{x} \\ \text{s. a. } & \boldsymbol{A}_{eq}\boldsymbol{x}=\boldsymbol{b}_{eq} \\ & \boldsymbol{A}\boldsymbol{x}\leq\boldsymbol{b}. \end{array} \end{equation}
Además, aprendimos a resolver los problemas en esta forma con la función linprog del paquete pyomo_utilities.py, proporcionando únicamente los parámetros $\boldsymbol{f}$, $\boldsymbol{A}$ y $\boldsymbol{b}$ ($\boldsymbol{A}_{eq}$ y $\boldsymbol{b}_{eq}$, de ser necesario).
In [16]:
f = -np.arange(1, 5)
A = np.array([[4, 3, 2, 1], [-1, -1, -1, -1]])
A = np.concatenate((A, -np.eye(4)), axis = 0)
b = np.array([10, -1, 0, 0, 0, 0])
Aeq = np.array([[1, 0, -1, 2]])
beq = np.array([2])
In [17]:
import pyomo_utilities
In [18]:
x, obj = pyomo_utilities.linprog(f, A, b, Aeq, beq)
In [19]:
x
Out[19]:
In [20]:
obj = -obj + 5
obj
Out[20]:
In [21]:
import pandas as pd
In [22]:
df = pd.DataFrame(columns=['Energia', 'Proteina', 'Calcio', 'Precio', 'Limite_diario'], index=['Avena', 'Pollo', 'Huevos', 'Leche', 'Pastel', 'Frijoles_cerdo'])
df.loc[:,'Energia']=[110, 205, 160, 160, 420, 260]
df.loc[:,'Proteina']=[4, 32, 13, 8, 4, 14]
df.loc[:,'Calcio']=[2, 12, 54, 285, 22, 80]
df.loc[:,'Precio']=[3, 24, 13, 9, 20, 19]
df.loc[:,'Limite_diario']=[4, 3, 2, 8, 2, 2]
df
Out[22]:
Luego de consultar expertos en nutrición tenemos que una dieta satisfactoria tiene por lo menos $2000$ kcal de energía, $55$ g de proteina, y $800$ mg de calcio.
Para imponer la variedad se ha decidido limitar el número de porciones diarias de cada una de las comidas como se indica en la tabla.
In [23]:
f = np.array(df.loc[:,'Precio'])
A = -np.array([df.loc[:,'Energia'], df.loc[:,'Proteina'], df.loc[:,'Calcio']])
A = np.concatenate((A, np.eye(6)), axis = 0)
A = np.concatenate((A, -np.eye(6)), axis = 0)
b = np.array([-2000, -55, -800])
b = np.concatenate((b, df.loc[:,'Limite_diario']))
b = np.concatenate((b, np.zeros((6,))))
In [24]:
x, obj = pyomo_utilities.linprog(f, A, b)
In [25]:
x
Out[25]:
In [26]:
obj
Out[26]:
El archivo forest_mex.csv contiene información histórica anual del porcentaje de área forestal de México. La primer columna corresponde a los años y la segunda corresponde al porcentaje de área forestal.
Tomado de: https://data.worldbank.org/indicator/AG.LND.FRST.ZS?view=chart.
Usando los años como variable independiente $x$ y el porcentaje de área forestal como variable dependiente $y$, ajustar polinomios de grado 1 hasta grado 3.
Mostrar en un solo gráfico los datos de porcentaje de área forestal contra los años, y los polinomios ajustados.
Graficar el error cuadrático acumulado contra el número de términos. ¿Cuál es el polinomio que mejor se ajusta?
Con los polinomios ajustados en el punto anterior, estime el año en que México quedará sin area forestal (suponiendo que todo sigue igual).
In [27]:
data_file = 'forest_mex.csv'
data = pd.read_csv(data_file, header = None)
In [28]:
x = data.iloc[:, 0].values
y = data.iloc[:, 1].values
plt.figure()
plt.plot(x, y)
plt.show()
In [29]:
beta1 = pyomo_utilities.curve_polyfit(x, y, 1)
beta2 = pyomo_utilities.curve_polyfit(x, y, 2)
beta3 = pyomo_utilities.curve_polyfit(x, y, 3)
In [30]:
yhat1 = beta1.dot(np.array([x**i for i in range(2)]))
yhat2 = beta2.dot(np.array([x**i for i in range(3)]))
yhat3 = beta3.dot(np.array([x**i for i in range(4)]))
plt.figure(figsize = (8,6))
plt.plot(x, y, '*b', label = 'datos reales')
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.show()
In [31]:
ems = []
ems.append(sum((y-yhat1)**2))
ems.append(sum((y-yhat2)**2))
ems.append(sum((y-yhat3)**2))
plt.figure(figsize = (8,6))
plt.plot(np.arange(3)+1, ems, '*b')
plt.xlabel('$n$')
plt.ylabel('$e_{ms}(n)$')
plt.show()
El polinomio que mejor se ajusta es el cúbico.
In [32]:
sym.var('a', real = True)
Af1 = beta1[0]+beta1[1]*a
Af2 = beta2[0]+beta2[1]*a+beta2[2]*a**2
Af3 = beta3[0]+beta3[1]*a+beta3[2]*a**2+beta3[3]*a**3
Af1, Af2, Af3
Out[32]:
In [33]:
a1 = sym.solve(Af1, a)
a2 = sym.solve(Af2, a)
a3 = sym.solve(Af3, a)
In [34]:
a1
Out[34]:
In [35]:
a2
Out[35]:
In [36]:
a3
Out[36]:
Hasta ahora hemos visto como diseñar clasificadores binarios dados un conjunto de entrenamiento. Sin embargo, no le hemos dado uso.
Después del diseño de un clasificador, lo que sigue es entregarle datos de entrada y que él los clasifique.
Para los datos de la tarea de clasificación binaria, diseñaremos un clasificador binario por regresión logística lineal utilizando únicamente los primeros 80 datos.
Luego, usaremos el clasificador diseñado para clasificar a los 20 datos restantes. ¿Cuántos datos se clasifican bien? ¿Cuántos mal?
In [37]:
X = 10*np.random.random((100, 2))
Y = (X[:, 1] > X[:, 0]**2)*1
In [38]:
plt.figure(figsize = (8,6))
plt.scatter(X[:,0], X[:,1], c=Y)
plt.show()
In [39]:
Xe = X[0:80, :]
Ye = Y[0:80]
In [40]:
B = pyomo_utilities.logreg_clas(Xe, Ye)
In [ ]:
In [ ]: