Ejemplo mínimo de aprendizaje automático¶
Métodos potenciales de prospección, FCAG, 2023.
En esta breve experiencia, empleamos TensorFlow. Para instalar el software en W*ndows con Miniconda instalado, escribir en la terminal de Miniconda:
conda> pip install --upgrade pip
conda> pip install tensorflow
A continuación, dado un conjunto $\{x_i,y_i\}$ de puntos de estación $x_i$ y observaciones $y_i$, más un modelo directo $f(x;a,b)$, intentaremos invertir el valor de los parámetros $a$ y $b$ que mejor ajustan las observaciones. La estructura del dato es similar a lo que veremos en las prácticas para un modelo directo de anomalÃa de gravedad.
$$\{x_i,\color{green}{y_i}\}, \quad f(x;{\color{blue}a},{\color{blue}b}) = e^{-{\color{blue}a}x^2} + \color{blue} b, \quad \color{green}{y_i}\approx f(x_i;\color{blue}a,\color{blue}b) = {\hat{\color{red}{y_i}}}.$$
import tensorflow as tf
import matplotlib.pyplot as plt
# Estilo para las notebooks preliminares:
plt.xkcd();plt.rcParams.update({'font.size': 16}) # tamaño de fuente
# El modelo y las observaciones
# A = 3.0
# B = 2.0
A0 = tf.constant(3.0)
B0 = tf.constant(2.0)
N = 100 # número de observaciones
# posiciones de las observaciones:
x = tf.random.normal(shape=[N])
# Ruido:
noise = 0.01*tf.random.normal(shape=[N])
# Calculamos y
#y = x * A + B + noise
#y = noise + tf.math.exp(-tf.math.pow(x,2) * A0) + B0
y = tf.math.exp(-tf.math.pow(x,2) * A0) + B0
# Visualizamos:
plt.figure(figsize=(7,5))
plt.scatter(x,y,label="dato $y_i$",color="royalblue",alpha=0.8)
plt.xlabel("x")
plt.ylabel("y")
plt.ylim([0.5,3.5])
plt.legend()
plt.show()
Modelo¶
Para realizar la inversión, haremos uso de tensorflow.fit(). Para ello escribimos el modelo directo como un objeto de la clase tensorflow.keras.Model, donde alcanza con detallar dos procesos. El proceso __init__() (con el cual inicializamos los parámetros del modelo directo) y el proceso call(), donde aplicamos el modelo al dato de entrada x.
class MyModel(tf.keras.Model):
def __init__(self, **kwargs):
super().__init__(**kwargs)
# Inicialización de los parámetros.
# En la práctica real, deberÃan ser inicializados aleatoriamente y de manera adecuada.
self.a = tf.Variable(initial_value = 0.0)
self.b = tf.Variable(initial_value = 0.0)
def call(self, x):
#return self.a * x + self.b
return tf.math.exp(-tf.math.pow(x,2) * self.a) + self.b
model = MyModel() # creamos una instancia del modelo
# Visualizamos:
plt.figure(figsize=(7,5))
plt.scatter(x, y, label="$y_i$", color="royalblue",alpha=0.8)
plt.scatter(x, model(x), label="$\hat y_i$ modelo inicial", color="tomato")
plt.xlabel("x")
plt.ylabel("y")
plt.ylim([0.5,3.5])
plt.legend()
plt.show()
Como observamos, el modelo inicial no resulta en un buen ajuste al dato. De manera cuantitativa, el error cuadrático medio (RMS) entre los datos y el resultado del modelo es:
def RMS(y,y0):
rms_2_per_sample = tf.square(y-y0) # RMS^2, por cada muestra individual de entrada
rms_2_sample_mean = tf.reduce_mean(rms_2_per_sample) # RMS^2, promediado en los datos de entrada
return tf.sqrt(rms_2_sample_mean).numpy() # RMS
print("RMS = ", RMS(model(x),y))
Entrenamiento¶
Definimos una función de costo (el cuadrado de la norma L2) y nos decidimos por un algoritmo optimizador.
- El optimizador es RMSprop.
- Utilizaremos un tamaño de paso o parámetro de aprendizaje (learning rate) fijo.
- Estimamos el valor del gradiente tomando el dato de entrada en grupos de
batch_sizemuestras (batch size). - Haremos un número de
epochspasadas (epochs) sobre el dato original completo.
Otros optimizadores a probar son descenso de gradiente estocástico (SGD) y Adam (adaptive moment estimation). El algoritmo Adam necesitará un número menor de iteraciones respecto a SGD para llegar a un resultado aceptable.
epochs = 500
batch_size = len(x)//10
learning_rate = 0.01
print("Épocas: ",epochs)
print("Tamaño de batch: ", batch_size)
print("Parámetro de aprendizaje: ", learning_rate)
# Función de costo:
@tf.function
def Myloss(target_y, predicted_y):
per_sample_loss = tf.square(target_y - predicted_y) # función de costo por cada par (y*_i,y_i)
return tf.reduce_mean(per_sample_loss) # función de costo promediada
# Optimizador:
#optimizer = tf.keras.optimizers.SGD(learning_rate = learning_rate)
#optimizer = tf.keras.optimizers.Adam(learning_rate = learning_rate)
optimizer = tf.keras.optimizers.RMSprop(learning_rate = learning_rate)
model = MyModel() # creamos una instancia del modelo
# `compile` define los parámetros de entrenamiento:
model.compile(optimizer = optimizer, loss = Myloss)
# Resumen del número de parámetros a invertir (`trainable parameters`):
model.build([len(x)]) # "build" para avisar la dimensión de la entrada asà vemos el summary y count_parameters
model.summary()
print("Número de parámetros del modelo: ", model.count_params())
# Entrenamiento
history = model.fit(x, y, epochs= epochs, batch_size= batch_size, verbose= 0) # para ver detalles al entrenar: verbose = 1
# Visualizamos
plt.figure(figsize=(7,5))
plt.scatter(x,y,label="$y_i$",color="royalblue",alpha=1)
plt.scatter(x,model(x),label="$\hat y_i$ modelo entrenado",color="orange",alpha=0.65)
plt.ylim([0.5,3.5])
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.show()
Apreciamos un mejor ajuste a los datos observados que el obtenido para el modelo inicial, con un RMS dado por:
print("RMS = ", RMS(model(x),y))
Los parámetros invertidos son:
a,b = model.weights
print("a* = ", a.numpy() )
print("b* = ", b.numpy() )
Por último, observamos la función de costo (normalizada en la primer época) en función del número de épocas.
# Visualizamos el costo (normalizado) a lo largo de las iteraciones:
perdida = history.history["loss"]
perdida = tf.math.divide(perdida,perdida[0])
plt.figure()
plt.plot(perdida,color="limegreen",lw=4)
#plt.plot(range(0,epochs,10),perdida[::10],"-o",color="limegreen",lw=0.5,ms=5)
plt.xlabel("# Épocas")
plt.ylabel("Costo (normalizado) []")
plt.show()
print("Costo actual: %1.6f" % Myloss(model(x), y).numpy())
Con estos conceptos y herramientas, en otra actividad práctica probaremos invertir datos de anomalÃa de gravedad para estimar parámetros fÃsicos y geométricos de un modelo directo no lineal.
Estimación de la ambiguedad¶
Corremos varias veces la inversión partiendo con distintos valores iniciales para los parámetros del modelo.
class MyModel(tf.keras.Model):
def __init__(self, **kwargs):
super().__init__(**kwargs)
# Inicialización de los parámetros.
# Aquà hacemos un poquito de trampa para no tener que invertir tantas épocas: generamos una
# distribución Gaussiana con valor medio en los valores reales de los parámetros a invertir
self.a = tf.Variable(tf.random.normal(shape=[1], mean = A0, stddev = 1.0, dtype=tf.dtypes.float32))
self.b = tf.Variable(tf.random.normal(shape=[1], mean = B0, stddev = 1.0, dtype=tf.dtypes.float32))
def call(self, x):
#return self.a * x + self.b
return tf.math.exp(-tf.math.pow(x,2) * self.a) + self.b
niter = 30
epochs = 100 # utilizamos un número pequeño de épocas para obtener resultados rápidos en esta demonstración.
results_list = []
for i in range(niter):
print(".",end =" ")
model = MyModel();
model.compile(optimizer = optimizer, loss = Myloss);
model.fit(x,y, epochs = epochs, batch_size = batch_size, verbose = 0) # verbose = 1 para ver detalles por época
results_list.append([model.a.numpy(), model.b.numpy()])
results = tf.stack(results_list)
# Valor medio y desviación estándar de cada parámetro:
As = results[:,0,:]
Bs = results[:,1,:]
Ainv_mean = tf.math.reduce_mean(As).numpy()
Ainv_std = tf.math.reduce_std(As).numpy()
Binv_mean = tf.reduce_mean(Bs).numpy()
Binv_std = tf.math.reduce_std(Bs).numpy()
# Gráficos
plt.figure()
plt.plot(As,":o",label="a*")
plt.plot(Bs,"-o",label="b*")
plt.xlabel("# n")
plt.ylabel("Parámetros invertidos")
plt.legend()
#plt.ylim([-tf.minimum(tf.abs(A),tf.abs(B)),+2*tf.maximum(tf.abs(A),tf.abs(B))])
plt.show()
# Valor medio y desviación estándar de cada parámetro:
print("a = {0}; <a*> = {1} +/- {2}".format(A0,Ainv_mean, Ainv_std))
print("b = {0}; <b*> = {1} +/- {2}".format(B0,Binv_mean, Binv_std))
Eso es todo por hoy.