# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Ajuste de curvas
#
# <img style="float: right; margin: 0px 0px 15px 15px;" src="https://upload.wikimedia.org/wikipedia/commons/a/a8/Regression_pic_assymetrique.gif" width="400px" height="125px" />
#
# > 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**:
# - https://en.wikipedia.org/wiki/Curve_fitting
# ___

# ## 0. Introducción
#
# Consideremos un polinomio de grado uno:
#
# $$y = \beta_1 x + \beta_0.$$
#
# Esta es una **línea recta** que tiene pendiente $\beta_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**
# 1. Encontrar la línea recta que pasa exactamente por los puntos $(0,1)$ y $(1,0)$.
# 2. Encontrar la parábola que pasa exactamente por los puntos $(-1,1)$, $(0,0)$ y $(1,1)$.
#
# **Solución**
# 1. Consideramos $y=\beta_1 x + \beta_0$. Evaluando en el punto $(0,1)$, obtenemos $\beta_1(0) + \beta_0 = 1$. Ahora, evaluando en el punto $(1,0)$, obtenemos $\beta_1(1) + \beta_0 = 0$. De esta manera,
# $$\left[\begin{array}{cc} 1 & 0 \\ 1 & 1\end{array}\right]\left[\begin{array}{c} \beta_0 \\ \beta_1\end{array}\right]=\left[\begin{array}{c} 1 \\ 0\end{array}\right].$$
# Resolviendo, $\beta_0=-\beta_1=1$.

# Importar numpy y el matplotlib.pyplot
import numpy as np
import matplotlib.pyplot as plt
# %matplotlib inline

# +
# Encontrar \beta_0 y \beta_1 resolviendo el sistema
A = np.array([[1, 0],
              [1, 1]])
b = np.array([1, 0])

beta = np.linalg.inv(A).dot(b)
#beta = np.dot(np.linalg.inv(A), b) 
beta
# -

# Graficar la recta encontrada junto con los puntos
x = np.array([0, 1])
y = np.array([1, 0])
plt.figure(figsize=(8,6))
plt.plot(x, y, 'ro', ms=5, label='Puntos')
plt.plot(x, beta[0]+beta[1]*x, label='Recta ajustada')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc='best')
plt.grid()

# 2. Consideramos $y=\beta_2 x^2 + \beta_1 x + \beta_0$. Evaluando en el punto $(-1,1)$, obtenemos $\beta_2(-1)^2 + \beta_1(-1) + \beta_0 = 1$. Ahora, evaluando en el punto $(0,0)$, obtenemos $\beta_2(0)^2 + \beta_1(0) + \beta_0 = 0$. Finalmente, evaluando en el punto $(1,1)$, obtenemos $\beta_2(1)^2 + \beta_1(1) + \beta_0 = 1$. De esta manera,
# $$\left[\begin{array}{ccc} 1 & -1 & 1 \\ 1 & 0 & 0 \\ 1 & 1 & 1 \end{array}\right]\left[\begin{array}{c} \beta_0 \\ \beta_1 \\ \beta_2 \end{array}\right]=\left[\begin{array}{c} 1 \\ 0 \\ 1 \end{array}\right].$$
# Resolviendo, $\beta_0=\beta_1=0$ y $\beta_2=1$.

# Encontrar \beta_0, \beta_1 y \beta_2
A = np.array([[1, -1, 1],
              [1, 0, 0],
              [1, 1, 1]])
b = np.array([1, 0, 1])

beta = np.linalg.inv(A).dot(b)
beta

# Graficar la parabola junto con los puntos
x = np.array([-1, 0, 1])
y = np.array([1, 0, 1])
xpuntos = np.linspace(-1,1,50)
plt.figure(figsize=(8,6))
plt.plot(x, y, 'ro', ms=5, label='Puntos')
plt.plot(xpuntos, beta[0]+beta[1]*xpuntos+beta[2]*xpuntos**2, label='Parabola ajustada')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc='best')
plt.grid()

# ### ¿Qué tienen en común los anteriores problemas?
# 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...

# Crear un conjunto de puntos ruidosos a partir de una recta
x = np.linspace(0, 1, 50)
y = 2+10*x+np.random.randn(50)

# Graficar
plt.figure(figsize=(8,6))
plt.plot(x, y, '*')
plt.xlabel('x')
plt.ylabel('y')
plt.grid()

# ### ¿Cómo ajustamos una curva a esto?

# ## 1. Problema básico
#
# <img style="float: right; margin: 0px 0px 15px 15px;" src="https://upload.wikimedia.org/wikipedia/commons/3/3a/Linear_regression.svg" width="400px" height="125px" />
#
# Consideramos que tenemos un conjunto de n pares ordenados de datos $(x_i,y_i)$, para $i=1,2,3,\dots,n$.
#
# ### ¿Cuál es la recta que mejor se ajusta a estos datos?
# 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]^T$, 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$$

# Notar que el problema anterior no es de programación lineal, ¿porqué?
#
# Para llevar a cabo la anterior minimización, la librería `SciPy` en su módulo `optimize` contiene la función `minimize`.

# Importar el módulo optimize de la librería scipy
import scipy.optimize as opt

# Función minimize
help(opt.minimize)


# Parámetros importantes:
# - fun: función $f(x)$, se debe definir antes de llamar minimize, como `def f(x): ... return ...`
# - x0: valor inicial. En una función no lineal, en general, hay múltiples mínimos. Dependiendo de la semilla caerá en uno de esos mínimos. Se ingresa como $x0 = \text{np.array}([x_{01},\dots,x_{0n}])$.
# - bounds: como en linprog.
# - constraints: funciones que definen las restricciones $g_i(x)$ y $h_j(x)$. Se definen igual que $f(x)$ y se ingresan como {'ineq': g_i, 'eq': h_j}.

# Primero debemos construir la función objetivo y la semilla inicial:

# +
# Definir funcion objetivo y punto inicial
def obj(beta, x, y):
    return np.sum((y-beta[0]-beta[1]*x)**2)

beta_inicial = [0, 0]
# -

resultado = opt.minimize(obj, beta_inicial, args=(x,y))

# Mostrar
resultado

# ¿Qué tan bien luce el ajuste?

# Coeficientes \beta_0 y \beta_1
beta = resultado.x
beta

# Grafica de los puntos y la recta ajustada
plt.figure(figsize=(8,6))
plt.plot(x, y, '*', label='Puntos')
plt.plot(x, beta[0]+beta[1]*x, label='Recta ajustada')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc='best')
plt.grid()

# Note que la pendiente es aproximadamente $10$ y el intercepto es aproximadamente $2$.
#
# La anterior idea se puede extender a ajuste polinomial...

# ## 2. Ajuste polinomial
#
# Ahora, considere el siguiente conjunto de datos...

# +
# Generamos 100 puntos ruidosos a partir de una senoidal
n = 100
x = np.linspace(np.pi/6, 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()


# -

# ### 2.1. ¿Se ajustará bien una recta?

# +
# Definir funcion objetivo y semilla
def obj(beta, x, y):
    return np.sum((y-beta[0]-beta[1]*x)**2)

beta_inicial = [0, 0]
# -

# Resolver
resultado = opt.minimize(obj, beta_inicial, args = (x,y))
resultado

# **Veamos $\beta$ para el ajuste con recta**

# Mostrar coeficientes
beta = resultado.x
beta

# Graficar
plt.figure(figsize=(8,6))
plt.plot(x, y, '*', label='Puntos')
plt.plot(x, beta[0]+beta[1]*x, label='Recta ajustada')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc='best')
plt.grid()


# ### 2.2. La recta no es buen ajuste... ¿Se ajustará bien una parabola?

# +
# Definir funcion objetivo y semilla
def obj2(beta, x, y):
    return np.sum((y-beta[0]-beta[1]*x-beta[2]*x**2)**2)

beta_inicial = [0, 0, 0]
# -

# Resolver
resultado = opt.minimize(obj2, beta_inicial, args=(x,y))
resultado

# **Veamos $\beta$ para el ajuste con parábola**

# Mostrar coeficientes
beta2 = resultado.x
beta2

# Graficar recta y parabola ajustadas
plt.figure(figsize=(8,6))
plt.plot(x, y, '*', label='Puntos')
plt.plot(x, beta[0]+beta[1]*x, label='Recta ajustada')
plt.plot(x, beta2[0]+beta2[1]*x+beta2[2]*x**2, label='Parabola ajustada')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc='best')
plt.grid()


# ### 2.3. Tampoco. Quizá un polinomio cúbico...

# +
# Definir funcion objetivo y semilla
def obj3(beta, x, y):
    return np.sum((y-beta[0]-beta[1]*x-beta[2]*x**2-beta[3]*x**3)**2)

beta_inicial = [0, 0, 0, 0]
# -

# Resolver
resultado = opt.minimize(obj3, beta_inicial, args=(x,y))
resultado

# **Veamos $\beta$ para el ajuste con cúbica**

# Mostrar coeficientes
beta3 = resultado.x
beta3

# Graficar recta, parabola y cubica
plt.figure(figsize=(8,6))
plt.plot(x, y, '*', label='Puntos')
plt.plot(x, beta[0]+beta[1]*x, label='Recta ajustada')
plt.plot(x, beta2[0]+beta2[1]*x+beta2[2]*x**2, label='Parabola ajustada')
plt.plot(x, beta3[0]+beta3[1]*x+beta3[2]*x**2+beta3[3]*x**3, label='Cubica ajustada')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc='best')
plt.grid()


# Mucho mejor. Entonces, ¿mientras más se suba el orden mejor la aproximación?
#
# ### 2.4. Ajustemos un polinomio de grado 7...

# +
# Definimos funcion objetivo y semilla
def obj7(beta, x, y):
    return np.sum((y-beta.dot([x**i for i in range(8)]))**2)

beta_inicial = np.zeros(8)
# -

# Resolvemos
resultado = opt.minimize(obj7, beta_inicial, args=(x,y))
resultado

# **De nuevo, veamos $\beta$**

# Mostrar coeficientes
beta7 = resultado.x
beta7

# **¡Cuidado! OVERFITTING...**
#
# Observar el tamaño de algunos coeficientes. Cuando los coeficientes son grandes, ¿qué pasa?

# Grafica de ajustes
plt.figure(figsize=(8,6))
plt.plot(x, y, '*', label='Puntos')
plt.plot(x, beta[0]+beta[1]*x, label='Recta ajustada')
plt.plot(x, beta2[0]+beta2[1]*x+beta2[2]*x**2, label='Parabola ajustada')
plt.plot(x, beta3[0]+beta3[1]*x+beta3[2]*x**2+beta3[3]*x**3, label='Cubica ajustada')
plt.plot(x, beta7.dot([x**i for i in range(8)]), label='Polinomio grado 7 ajustada')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc='best')
plt.grid()

# Es conveniente ver el error como función del orden del polinomio... **selección de modelos**

# +
# Error cuadratico
e_ms = []
def obj(b, x, y, n):
    return np.sum((y-b.dot([x**i for i in range(n+1)]))**2)

for i in range(7):
    b0 = np.zeros((i+2,))
    res = opt.minimize(obj, b0, args=(x,y,i+1))
    e_ms.append(res.fun)
    
plt.figure(figsize=(8,5))
plt.plot(np.arange(7)+1, e_ms, 'o')
plt.xlabel('orden')
plt.ylabel('error')
plt.show()


# -

# En efecto, parece que con $3$ es suficiente.

# ### ¿Cómo prevenir el *overfitting* sin importar el orden del modelo?

# ## 3. Regularización
#
# 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. 

# ### 3.1. Ridge
#
# $$\boldsymbol{\beta}^{ridge} = \arg \min_{\boldsymbol{\beta}} \left|\left|\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right|\right|^2 + \lambda\left|\left|\boldsymbol{\beta}\right|\right|^2$$

# +
def obj_ridge(b, x, y, n, l):
    return np.sum((y-b.dot([x**i for i in range(n+1)]))**2)+l*np.linalg.norm(b)**2

b0 = np.random.random((8,))
res = opt.minimize(obj_ridge, b0, args=(x,y,7,0.1))
yhat7_ridge = np.array([x**j for j in range(8)]).T.dot(res.x)

plt.figure(figsize=(8,6))
plt.plot(x, y, '*b', label = 'datos')
plt.plot(x, beta7.dot([x**i for i in range(8)]), '-c', label = 'ajuste 7')
plt.plot(x, yhat7_ridge, '--r', label = 'ajuste 7_ridge')
plt.legend(loc = 'best')
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.grid()
plt.show()
# -

res.x


# ### 3.2. Lasso
#
# $$\boldsymbol{\beta}^{lasso} = \arg \min_{\boldsymbol{\beta}} \left|\left|\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right|\right|^2 + \lambda\left|\left|\boldsymbol{\beta}\right|\right|_1$$
#
# 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|$.

# +
def obj_lasso(b, x, y, n, l):
    return np.sum((y-b.dot([x**i for i in range(n+1)]))**2)+l*np.linalg.norm(b,1)

b0 = np.random.random((8,))
res = opt.minimize(obj_lasso, b0, args=(x,y,7,0.1))
yhat7_lasso = np.array([x**j for j in range(8)]).T.dot(res.x)

plt.figure(figsize=(8,6))
plt.plot(x, y, '*b', label = 'datos')
plt.plot(x, beta7.dot([x**i for i in range(8)]), '-c', label = 'ajuste 7')
plt.plot(x, yhat7_ridge, '--r', label = 'ajuste 7_ridge')
plt.plot(x, yhat7_lasso, '--k', label = 'ajuste 7_lasso')
plt.legend(loc = 'best')
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.grid()
plt.show()
# -

res.x

# ## 4. Ajuste robusto
#
# Ahora, consideremos de nuevo el caso de la línea recta con un par de puntos atípicos al inicio y al final...

# +
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()
# -

# Solucionamos el problema normalmente...

# +
b0 = np.random.random((2,))
res = opt.minimize(obj, b0, args=(x,y,1))
yhat = np.array([x**j for j in range(2)]).T.dot(res.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()
# -

res.x


# 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*.

# +
def huber(a, d):
    if np.abs(a)<=d:
        return a**2
    else:
        return d*(2*np.abs(a)-d)

def obj_robust(b, x, y, n, d):
    return np.sum(np.vectorize(huber)(y-b.dot([x**i for i in range(n+1)]), 1.345))

b0 = np.random.random((2,))
res = opt.minimize(obj_robust, b0, args=(x,y,1,1.345))
yhat = np.array([x**j for j in range(2)]).T.dot(res.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()
# -

res.x


# Mejor...

# ## 5. Actividad
#
# 1. Ajustar polinomios de grado 1 hasta grado 7 a los siguientes datos.
# 2. Graficar el error cuadrático acumulado contra el número de términos, y elegir un polinomio que ajuste bien y su grado no sea muy alto.
# 3. Para el grado de polinomio elegido, realizar el ajuste con ridge con coeficiente de 0.01.
# 4. Comparar los beta.
#
# Abrir un nuevo notebook, llamado `Tarea5_ApellidoNombre` y subirlo a moodle en el espacio habilitado. Tarea para el Lunes 17 a las 23:00.

def f(x):
    return np.exp(-x**2/2)/np.sqrt(2*np.pi)


# +
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()
# -

# # Avisos:
#
# ## Recordar tarea para hoy y tarea para el jueves.
#
# ## Evaluación primer módulo: (Por definir)
# ### Se las entrego el Lunes 17 de Septiembre (dentro de una semana) por moodle y tienen hasta el Viernes 21 de Septiembre.
#
# ## Proyecto:
# ### 1. Elegir integrantes para proyecto. Mínimo 2, máximo 3 (sin excepción). Entregarme una hoja por equipo con los nombres de los integrantes ya.
# ### 2. Deben elegir un tema para proyecto que se pueda resolver como un problema de optimización (preferiblemente, relacionado con su carrera). 
# ### 3. Para la siguiente semana, a más tardar, deben acercarse a mi con su tema de proyecto. Juntos, definiremos el alncance.
# ### 4. Fecha de entrega y presentación: Lunes 24 de septiembre.

# <script>
#   $(document).ready(function(){
#     $('div.prompt').hide();
#     $('div.back-to-top').hide();
#     $('nav#menubar').hide();
#     $('.breadcrumb').hide();
#     $('.hidden-print').hide();
#   });
# </script>
#
# <footer id="attribution" style="float:right; color:#808080; background:#fff;">
# Created with Jupyter by <NAME>.
# </footer>

StarcoderdataJupyterScriptsDedupFiltered