In [ ]:
"""
IPython Notebook v4.0 para python 2.7
Librerías adicionales: numpy, matplotlib
Contenido bajo licencia CC-BY 4.0. Código bajo licencia MIT. (c) Sebastian Flores.
"""
# Configuracion para recargar módulos y librerías
%reload_ext autoreload
%autoreload 2
%matplotlib inline
from IPython.core.display import HTML
import warnings
warnings.filterwarnings("ignore")
HTML(open("style/mat281.css", "r").read())
Antes: Modelo de publicación actual privilegia resultados exitosos, por sobre documentación de experimentos fallidos. Existe un sesgo de publicación (publication bias).
Durante: Nunca hay tiempo suficiente. Priorización del tiempo a obtener más y mejores resultados, no a examinar casos fallidos.
Después: Posible y fácil, pero rara vez realizado, puesto que es opcional y requiere fuerza de voluntad. Común en la milicia, bajo el nombre de debriefing.
Certámenes y Tareas:
Proyectos en grupo:
Debriefing no altera resultados pasados, pero permite mejorar incrementalmente y detectar errores metodológicos.
Las mediciones son siempre indirectas. Lo que se mide no es necesariamente lo que se necesita.
Longitud: comparación directa (regla), comparación indirecta (sensor laser).
Temperatura: mediante dilatación (termómetro tradicional), cambio de densidad ( termómetro de Galileo), o cambio de resistencia eléctrica (termómetro digital).
Los parámetros rara vez se conocen con buena precisión, incluso cuando el modelo físico-matemático ha sido bien definido. Por ello, siempre es necesario tomar un enfoque conservador en la incerteza de los parámetros involucrados dependiendo de la aplicación del modelo.
Por ejemplo, en el caso de minería subterránea:
Si el problema involucra seguridad en los túneles mineros, se utiliza un modelo conservador que subestima los parámetros físicos de la roca: : la roca del modelo es menos resistente que en la realidad. Si el modelo indica que el túnel es resistente, tenemos confianza que en la realidad el tunel es seguro.
Si el problema involucra factibilidad de la técnica extractiva, se utiliza un modelo que sobreestima los parámetros físicos de la roca: la roca del modelo es más resistente que en la realidad. Si el modelo indica que la roca se quiebra adecuadamente durante el caving, tenemos confianza que la roca se quiebra durante el caving.
El análisis de los errores en ingeniería no se realiza etapa por etapa. Tradicionalmente y por conveniencia se supone que el error es únicamente experimental, es decir, se han detectado y corregido los errores ilegítimos o humanos.
Errores experimentales se dividen en 2 categorías:
El modelamiento no es un objetivo en sí mismo.
"All models are wrong,
some are useful" - George Box"Since all models are wrong the scientist
cannot obtain a correct one by excessive
elaboration. On the contrary following William
of Occam he should seek an economical
description of natural phenomena. [...]
the ability to devise simple but evocative
models is the signature
of the great scientist [...]"- George BoxSupongamos que ciertos datos provienen de la siguiente caja negra: $$ y(x) = 5 \cos\Big( \frac{\pi x}{2}\Big) + \mathcal{N}(0,1)$$
Tenemos 200 mediciones, y deseamos ajustar un modelo de tipo polinomial a los datos. Después de todo, sabemos que
$$ cos(x) = 1 - \frac{1}{2!}x^2 + \frac{1}{4!}x^4 - \frac{1}{6!}x^6 + \frac{1}{8!}x^8 + ...$$Además, 200 datos son bastantes, por lo que, ¿qué tan mal lo podríamos hacer?
In [ ]:
import numpy as np
from matplotlib import pyplot as plt
# Load data
x, y = np.loadtxt("data/data.txt").T
x_bb = np.sort(x)
y_bb = 5 * np.cos(0.5*np.pi*x_bb/2)
# Plot data
plt.figure(figsize=(16,8))
plt.plot(x, y, 'bs', alpha=0.5, label=u"Medición")
plt.plot(x_bb, y_bb, 'k-', lw=2.0, label=u"Relación determinista")
plt.xlim([-2.5,2.5])
plt.ylim([-5,10])
plt.xlabel("x []", fontsize=20)
plt.ylabel("y []", fontsize=20)
plt.legend(numpoints=1, loc="lower center", fontsize=20)
plt.show()
In [ ]:
# Define error function
def error(vector_e):
return abs(vector_e).mean()
# Parameters
nmax = 70 # max polynomial
# Load data
x, y = np.loadtxt("data/data.txt").T
# Fit best several models and record training error and prediction error
order_polynomial = range(1, nmax+1)
error_polynomial = []
for n in order_polynomial:
fit_n = np.polyfit(x, y, n) # Obtains the best fitted polynomial of degree n
pol_n = np.poly1d(fit_n) # Creates the polynomial with coefficients as in fit n
error_polynomial.append( error(y - pol_n(x)) )
# Plot data
plt.figure(figsize=(16,8))
plt.semilogy(order_polynomial, error_polynomial, "-sb", lw=2.0, label="Error")
plt.xlabel("x []", fontsize=20)
plt.ylabel("y []", fontsize=20)
plt.xlabel("Grado del polinomio", fontsize=20)
plt.ylabel("Error", fontsize=20)
plt.show()
In [ ]:
# Parameters
n = 1 # order of polynomial
# Load data
x, y = np.loadtxt("data/data.txt").T
# Fit best several models and record training error and prediction error
fit_n = np.polyfit(x, y, n) # Obtains the best fitted polynomial of degree n
pol_n = np.poly1d(fit_n) # Creates the polynomial with coefficients as in fit n
# Plot data
x_bb = np.sort(x)
y_bb = 5 * np.cos(0.5*np.pi*x_bb/2)
plt.figure(figsize=(16,8))
plt.plot(x, y, 'bs', alpha=0.5, label=u"Medición")
plt.plot(x_bb, y_bb, 'k-', lw=2.0, label=u"Relación determinista")
plt.plot(x_bb, pol_n(x_bb), 'r-', lw=2.0, label=u"Polinomio grado %d" %n)
plt.xlabel("x []", fontsize=20)
plt.ylabel("y []", fontsize=20)
plt.legend(numpoints=1, loc="lower center", fontsize=20)
plt.show()
¿Qué hacer? ¿Existe algún método que me permita seleccionar cuantitativamente el mejor orden?
¿Qué tal si entrenamos el modelo con algunos datos, y después comprobamos que tal se comporta el modelo en los datos reservados? (holdout set technique).
In [ ]:
# Parameters
n = 1 # order of polynomial
# Load data
x, y = np.loadtxt("data/data.txt").T
# Split into training and prediction data
t = int(len(x)*.7)
x_t = x[:t]
y_t = y[:t]
x_p = x[t:]
y_p = y[t:]
# Fit best several models and record training error and prediction error
fit_n = np.polyfit(x_t, y_t, n) # Obtains the best fitted polynomial of degree n
pol_n = np.poly1d(fit_n) # Creates the polynomial with coefficients as in fit n
# Plot data
x_bb = np.sort(x)
y_bb = 5 * np.cos(0.5*np.pi*x_bb/2)
plt.figure(figsize=(16,8))
plt.plot(x_t, y_t, 'bs', alpha=0.5, label=u"Training Set")
plt.plot(x_p, y_p, 'gs', alpha=0.5, label=u"Testing Set")
plt.plot(x_bb, y_bb, 'k-', lw=2.0, label=u"Relación determinista")
plt.plot(x_bb, pol_n(x_bb), 'r-', lw=2.0, label=u"Polinomio grado %d" %n)
plt.xlabel("x []", fontsize=20)
plt.ylabel("y []", fontsize=20)
plt.legend(numpoints=1, loc="lower center", fontsize=20)
plt.show()
In [ ]:
# Parameters
nmax = 70 # max polynomial
# Load data
x, y = np.loadtxt("data/data.txt").T
# Fit best several models and record training error and prediction error
order_polynomial = range(1, nmax+1)
error_t = []
error_p = []
for n in order_polynomial:
fit_n = np.polyfit(x_t, y_t, n) # Obtains the best fitted polynomial of degree n
pol_n = np.poly1d(fit_n) # Creates the polynomial with coefficients as in fit n
error_t.append( error(y_t - pol_n(x_t)) )
error_p.append( error(y_p - pol_n(x_p)) )
# Plot the errors
plt.figure(figsize=(16,8))
plt.semilogy(order_polynomial, error_t, "-sb", lw=2.0, label="Training error")
plt.semilogy(order_polynomial, error_p, "-og", lw=2.0, label="Prediction error")
plt.legend(numpoints=1, loc="lower center", fontsize=20)
plt.xlabel("Grado del polinomio", fontsize=20)
plt.ylabel("Error", fontsize=20)
plt.show()