Jorge III Altamirano-Astorga, Luz Aurora Hernández-Martínez, Ita-Andehui Santiago-Castillejos.

Profesor: Dr. Edgar Francisco Román-Rangel

Resumen

Desarrollaremos un proyecto de investigación basados en un sensor de la calidad del aire que tenemos dentro de casa de uno de los participantes con el fin de estudiar, analizar, explorar y entender su relación e influencia con los fenómenos externos (calidad del aire de la ciudad y variables atmosféricas) para poder predecir la calidad del aire en el interior de casa con las mediciones de compuestos orgánicos volátiles, los cuales tienen alta probabilidad de ser perjudiciales para la salud.

Introducción

El ITAM y muchos de sus miembros tenemos residencia en la Ciudad de México. Esta metrópolis es una de las más contaminadas en el continente y en el mundo. A causa de la pandemia, una gran cantidad de la población pasamos mucho de nuestro tiempo confinados en espacios cerrados, típicamente nuestras viviendas. Queremos saber cómo influyen los factores atmosféricos y la contaminación de la Ciudad en la contaminación de un espacio cerrado.

Hemos almacenado los datos de este sensor desde Febrero 2021 tratando de mantenerlos en un área común que no tiene ventilación directa para evitar perturbaciones en las lecturas y que sea influído directamente por la contaminación exterior, así como de otras fuentes de emisión (cocina).

Es importante destacar que este sensor no detecta contaminantes primarios, tales como: óxidos de nitrógeno (NOx), dióxido de carbono (CO2), monóxido de carbono (CO), ozono (O3), más bien mide los compuestos orgánicos volátiles, conocidos por el acrónimo anglosajón VOCs. Los VOCs típicamente son muchos de los olores que percibimos, los cuales son disoluciones de compuestos en el aire.

Estos compuestos orgánicos volátiles se han comprobado como nocivos a la salud y posibles cancerígenos, lo cual nos despertó el interés. Ejemplos de estos compuestos orgánicos son: el humo del cigarro, humo causado por cocinar alimentos, la utilización de agentes de limpieza (cloro y basados en amoniaco) y fuentes volátiles varias (como solventes, pinturas, quitaesmaltes), entre otros.

Fuentes de Datos y Variables

Tenemos los siguientes de fuentes de datos con las siguientes variables:

1. Sensor de Contaminantes en Interior de Casa ME680: contamos aproximadamente 2.1 millones de registros. Las lecturas del sensor son cada 3 segundos.

   • Temperatura: variable numérica en grados Celsius (C) con una resolución de $0.01C$ y una precisión de $\pm 0.5C$.

   • Presión: variable numérica en hectopascales (hPa) con una resolución de $0.18~hPa$ y una precisión de $\pm 0.12~hPA$.

   • Humedad: variable numérica en porcentaje de humedad relativa (%rH) con una resolución de $0.008~\%rH$ y una precisión de $\pm 3\%rH$.

   • Resistencia del Gas: variable numérica de la resistencia eléctrica opuesta al elemento sensible del sensor medida en Ohms.

   • IAQ: variable numérica medida en el índice de calidad del aire americano en interior (IAQI, aunque utilizaremos la nomenclatura IAQ) con una resolución de $1~ IAQ$. La precisión del sensor variable que no excede 5% se guarda en una variable independiente.

   • Precisión del sensor: variable categórica ordinal con valores en el rango de [0,3]:

     • 0: periodo de estabilización o no operativo.

     • 1-2: periodo operativo.

     • 3: precisión máxima y operación óptima.

   • Fecha y hora: variable numérica basado en UNIX/POSIX epoch que denota el tiempo desde el 01/01/1970 00:00:00.0 UTC. El tiempo está sincronizado por NTP al Centro Nacional de Metrología de México (Hora Oficial del País).

2. Datos Abiertos del Gobierno (SINAICA) la Calidad del Aire del Gobierno de la Ciudad de México con 2,170 observaciones. Los valores reportados son cada Hora, todos son variables continuas.

   • CO: Monóxido de Carbono medido en partes por millón (ppm).

   • NO, NO2, NOx: Familia de óxidos de nitrógeno. 3 variables correspondientes a la familia de óxidos de nitrógeno. Medidos en partes por billón (ppb).

   • O3: Ozono medido en partes por billón (ppb).

   • PM2.5, PM10: partículas microscópicas de 2.5 micras y 10 micras. Medido en microgramos por metro cúbico ($\mu g / m^3$)

   • SO2: dióxido de azufre. Medido en partes por billón (ppb).

   • Fecha y hora: entregado en zona horaria del Centro de México.

     Tomamos los datos de las estaciones meteorológicas cercanas en un radio de 10 km del sensor interior.

Todos los datos que tenemos fueron divididos en un set de entrenamiento y pruebas tanto para el desarrollo de los modelos predictivos, como para las imputaciones. Los datos de SINAICA se dividieron en conjuntos correspondientes a 70% de los datos y 30% para pruebas. Los datos del sensor, al ser más, utilizamos una división de 80-20.

Problemáticas

• Datos del Gobierno de la Ciudad de México: estos datos no están siendo actualizados de manera cotidiana, por lo que en ocasiones tienen atrasos en publicar la información actual. Ejemplo, al 31 de marzo no habían subido actualizaciones desde el 28 de enero. Se pudo contrarrestar este problema utilizando la fuente de datos federal (SINAICA) en esta sección: proyectofinal01_imputaciones_sinaica.html.

• No pudimos utilizar datos meteorológicos de terceros dado que nos son fuentes abiertas y tenían un costo los datos históricos.

• Precisión y manipulación de los datos del gobierno: observamos demasiadas observaciones faltantes en todas las estaciones meteorológicas, ya sea por manipulación directa para ocultar información o bien por descuido y falta de interés del gobierno para publicar estos datos.

• Estabilidad y precisión de la toma de registros en el sensor: tuvimos típicamente interrupciones menores a 1 minuto causado por actualizaciones de software y fallas eléctricas.

• Propiedad intelectual: el algoritmo del sensor es cerrado para convertir de la variable gasResistance a la variable IAQ. Sólo proporcionan un objeto binario, el cual procesa la variable gasResistance. [1] [2]

  • Al parecer esta variable se genera como una serie de tiempo y llega a guardar hasta 24 hrs de estados previos de ejecución del sensor. [3]

  • Este código es protegido como propiedad intelectual del fabricante del sensor (Bosch), lo cual impide que se pueda hacer desensamble o ingeniería inversa.

  • Se pudiera utilizar el objeto binario en un ejecutable en Linux o Windows, donde pudiese ser convertida no desde el sensor, sino desde nuestras predicciones de la variable gasResistance.

Este problema lo logramos superar robusteciendo la red neuronal para que tuviéramos las predicciones en la unidad IAQ que se describió en la sección anterior: Variables.

Análisis Exploratorio Inicial

Se realizó un análisis exploratorio de los datos y recopilación de los datos de fuentes externas. Estas se pueden encontrar en proyectofinal00_eda.html

A continuación mostramos algunas gráficas de las lecturas del sensor.

%matplotlib inline
from IPython.display import (display, Markdown, Image, 
                             clear_output, Latex, Math)
import pandas as pd
import re, os, sys, shelve, time, dill, io
from pickle import PicklingError
from dill import Pickler, Unpickler
shelve.Pickler = Pickler
shelve.Unpickler = Unpickler
import numpy as np
import matplotlib.pyplot as plt
from plotnine import *
import plotnine.options as p9opts
#figure_size = (6.4, 4.8)
p9opts.figure_size = (2.5, 1.2)
#p9opts.dpi = 84
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import MinMaxScaler

airdata = pd.read_pickle("data/airdata/air-imputated.pickle.gz")
display(Markdown(f"* Rango de fechas obtenidas: {airdata.datetime.min().strftime('%Y-%m-%d %H:%M:%S')} al {airdata.datetime.max().strftime('%Y-%m-%d %H:%M:%S')}"))
display(Markdown(f"* Número de registros: {airdata.shape[0]:3,}"))
display(Markdown(f"* Promedio de IAQ: {airdata.IAQ.mean():.2f} desviación estándar: {airdata.IAQ.std():.2f}"))

• Rango de fechas obtenidas: 2021-02-12 06:05:35 al 2021-04-24 22:16:20

• Número de registros: 2,068,354

• Promedio de IAQ: 161.23 desviación estándar: 72.85

_ = (
    ggplot(airdata.sample(frac=0.3, random_state=175904), 
           aes(x = "datetime", y = "temperature", color="pressure")) +
    geom_jitter(alpha=0.05) +
    theme(axis_text_x=element_text(angle=45), 
          plot_title=element_text(size=11)) +
    labs(x="Fecha", y="Temp (C)", color="Presión (hPa)", 
         title="Gráfica de Temperatura y\nPresión a lo Largo del Tiempo.") 
).draw()

#_ = (
#    ggplot(airdata.sample(frac=0.3, random_state=175904), 
#           aes(x = "datetime", y = "IAQ", color="iaqAccuracy")) +
#    geom_jitter(alpha=0.05, size=1.25) +
#    theme(axis_text_x=element_text(angle=45)) +
#    labs(x="Fecha", color="Precisión del Sensor",
#         title="Gráfica de IAQ y Precisión\ndel Sensor a lo Largo del Tiempo")
#).draw()

Modelos

Proponemos utilizar un modelo de aprendizaje profundo con distintos tipos de neuronas artificiales:

• DNN: Las redes neuronales densas son redes neuronales profundas que son la base para las redes neuronales artificiales (ANN) con múltiples capas entre las capas de entrada y salida.

• RNN: Las redes neuronales recurrentes son la arquitectura más clásica y utilizada para problemas de predicción de series temporales;

• CNN: red neuronal convolucional es muy popular en procesamiento de imágens, sin embargo es útil para modelos unidimensionales y en más de 2 dimensiones.

• LSTM: Las redes neuronales de "Memoria Largo-Corto Plazo" (LSTM) que son una evolución de las RNN desarrolladas para superar el problema del gradiente que desaparece;

• Mezcla de los mejores modelos, en este caso CNN + LSTM.

Utilizamos series de tiempo basándonos en el modelo de ventanas de tiempo donde se crearon matrices con tiempos pasados [4] [5]

Además de utilizar técnicas de series de tiempo, estadística frecuentista y bayesiana para el análisis de los datos y las imputaciones que fueron necesarias para tener esta secuencia de valores para la serie de tiempo. Exploramos las siguientes técnicas.

• KNN: K-Vecinos más cercanos.

• GLM: Métodos Lineales Generalizados: tanto frecuentista, como bayesiano con el método del muestreador de Gibbs.

• Medias

• Hot-Deck

• Interpolación con ruido.

Esta última fue la técnica que utilizamos por adecuarse mejor a la secuencia de los datos. Como se muestra en un ejemplo a continuación:

_ = airdata
_.index = airdata.datetime
_ = _["2021-02-18 23:10":"2021-02-18 23:35"].copy()
_["datetime2"] = _["datetime"].dt.strftime('%H:%M:%S')
#display(_)
_ = (
  ggplot(_) +
  geom_point(aes(x="datetime", y="humidity", color="imputated")) +
  scale_x_datetime(date_breaks='5 minute', date_labels="%H:%M:%S") +
  theme(axis_text_x=element_text(angle=45), 
          plot_title=element_text(size=11)) +
  labs(title="Gráfica Ejemplo de Datos Imputados.")
).draw()

Se pueden encontrar el detalle en la sección: proyectofinal01_imputaciones_airdata.html

Trabajos relacionados

Existen los siguientes trabajos relacionados:

1. Examen final para la materia de "Modelos de Gran Escala" con la Prof. Liliana Millán. Utilizamos los datos de las estaciones de monitoreo ambiental de la Ciudad de México y los datos de afluencia de las estaciones Ecobici. Buscamos establecer la relación entre la afluencia de las estaciones de Ecobici y la disminución (o aumento) de contaminantes en las inmediaciones a estas estaciones. Encontramos que las estaciones con alta afluencia tienen mayor contaminación, pero pudiera ser provocado porque la afluencia es para todos los medios de transporte, incluídos los emisores de contaminación. Fueron utilizados metodos de aprendizaje máquina.

2. Development of indoor environmental index: Air quality index and thermal comfort index [6] Es un estudio de la Universidad Tecnológica de Malasia. Utilizaron los datos de un conjunto de sensores, similares al nuestro, para establecer la relación entre contaminantes en interior (IAQ) y el índice de comodidad térmica (TCI) mediante un nuevo índice de calidad y comodidad ambiental de interiores mediante un monitoreo en "tiempo real". El modelo propuesto está basado en una suma ponderada.

3. Tensorflow Tutorial on Time-Series Forecasting [4]. Los creadores de Tensorflow, a manera de demostración, utilizaron los datos de Biogeoquímica del Instituto Max-Planck de Alemania para predecir el clima a partir de datos históricos utilizando redes neuronales en combinación con series de tiempo. Utilizaron redes neuronales de los tipos RNN y LSTM.

Solución

Preprocesamiento

Hicimos un preprocesamiento de los datos basado en los siguientes puntos:

• Conjunto de datos de Entrenamiento y Pruebas que fue descrito en la sección de Fuentes de Datos y Variables.

• Imputación: imputamos las observaciones faltantes con datos de la calidad del aire de la Ciudad de México, en específico con la información del monitor de la estación Camarones, esto considerando su distancia a la ubicación de nuestro sensor. Para hacer la imputación usamos interpolación, aunque también exploramos K Vecinos más Cercano y Métodos Lineales Generalizados, entre otros. También es importante mencionar que procesamos los datos como una serie de tiempo usando un tutorial de tensorflow y Keras [5].

Datos del Gobierno y del Sensor

En esta sección se puede encontrar el procesamiento de los datos de las estaciones de monitoreo de contaminantes de la Ciudad de México y los Datos del Sensor: proyectofinal03_sensor-sinaica.html

Modelo Propuesto

El modelo que mejores resultados nos dio fue un modelo de redes neuronales que combina en una sola red: subredes convolucionales (CNN), memoria largo-corto plazo (LSTM) y redes densas (DNN); observamos las fortalezas de cada una de dichas subredes, pero no observamos en gran medida sus debilidades. El modelo aprendió una representación interna de los datos de la serie temporal y logró el mejor rendimiento.

En el siguiente documento encontrará la arquitectura del modelo: proyectofinal04_mejora_modelo.html.

Resultados

A continuación presentaremos los desempeños desescalados con el fin de comparar el comportamiento de los modelos y poder decir la influencia y caapacidad de predicción:

try:
  from google.colab import drive
  drive.mount('/content/drive', force_remount=False)
  from google.colab import files
except:
  ;

base_url = ""
# File Loaders
try:
    base_url = "drive/MyDrive/Colab Notebooks/Final"
    uploaded = os.path.join(base_url, "data/sinaica-imputated.pickle.gz")
    if(not os.path.isfile(uploaded)):
        from google.colab import files
        uploaded = files.upload()
except:
    base_url = ""
    uploaded = "data/sinaica/sinaica-imputated.pickle.gz"

def render_mpl_table(data, col_width=3.0, row_height=0.625, font_size=14,
                     header_color='#40466e', row_colors=['#f1f1f2', 'w'], edge_color='w',
                     bbox=[0, 0, 1, 1], header_columns=0,
                     ax=None, **kwargs):
    """
    Taken from https://stackoverflow.com/a/39358722/7323086
    """
    if ax is None:
        size = (np.array(data.shape[::-1]) + np.array([0, 1])) * np.array([col_width, row_height])
        fig, ax = plt.subplots(figsize=size)
        ax.axis('off')

    mpl_table = ax.table(cellText=data.values, bbox=bbox, colLabels=data.columns, **kwargs)

    mpl_table.auto_set_font_size(False)
    mpl_table.set_fontsize(font_size)

    for k, cell in  six.iteritems(mpl_table._cells):
        cell.set_edgecolor(edge_color)
        if k[0] == 0 or k[1] < header_columns:
            cell.set_text_props(weight='bold', color='w', size=font_size*1.05)
            cell.set_facecolor(header_color)
        else:
            cell.set_facecolor(row_colors[k[0]%len(row_colors) ])
    plt.show()

#df.dropna(inplace=True)
clear_output()

def performance_plot(history, a=None, b=None, 
                    metrics=["accuracy", "val_accuracy"],
                    plot_validation=True,
                    title="Gráficas de Desempeño."):
  """
  Prints performance plot from a, to b on a history dict.
  
  Inputs:
  history: dict containing "loss" and "accuracy" keys
  a: epoch start
  b. last epoch
  metrics: plot these metrics (train and validation). Always 2.
  plot_validation: boolean indicating if validation data should be plotted.
  a: from this epoch
  b: to this epoch    
  """
  if a is None:
      a = 0
  if b is None:
      b = len(history['loss'])
  a = np.min((a,b))
  b = np.max((a,b))

  imgrows = (len(metrics) + 1) / 2
  imgrows = np.round(imgrows, 0)
  imgrows = int(imgrows)
  #print(imgrows)

  # Plot loss
  plt.figure(figsize=(14, 5
                      *imgrows))
  plt.suptitle(title)
  plt.subplot(imgrows, 2, 1)
  plt.title('Loss')
  plt.plot(history['loss'][a:b], label='Training', linewidth=2)
  if plot_validation:
    plt.plot(history['val_loss'][a:b], label='Validation', linewidth=2)
  plt.legend()
  plt.xlabel('Epoch')
  plt.ylabel(f'Loss')
  quantiles = np.quantile(range(a, b), 
                          [.2, .4, .6, .8]).round(0).astype(int)
  quantiles = np.insert(quantiles, 0, [a])
  quantiles += 1
  quantiles = np.append(quantiles, [b-1])
  plt.xticks(ticks=quantiles-a,
              labels=quantiles)
  plt.grid(True)

  # Plot accuracy
  for i, metric in enumerate(metrics): 
    #print(f"metric: {metric}, i: {i}")
    #print(f"mean metric: {np.mean(history[metric])}")
    plt.subplot(imgrows, 2, i+2)
    plt.title(metric)
    plt.plot(history[metric][a:b], label='Training', 
              linewidth=2)
    if plot_validation:
      plt.plot(history["val_" + metric][a:b], 
                label='Validation', linewidth=2)
    plt.legend()
    plt.xlabel('Epoch')
    plt.ylabel(metric)
    #plt.xlim(a, b)
    #print(range(0, b-a))
    plt.xticks(ticks=quantiles-a, 
                labels=quantiles)
    plt.grid(True)

  plt.show()

#render_mpl_table(df.head().applymap(shorten), col_width=5)

sinaica = pd.read_pickle(uploaded)
airdata = pd.read_pickle(os.path.join(base_url, "data/airdata/air-imputated.pickle.gz"))
#sinaica.head()

• Tabla Comparativa de los Datos del Gobierno y del Sensor

models = []
object_names = []
models_path = os.path.join(base_url, "models-sinaica")

for y in [x for x in os.listdir(models_path) if x.endswith("dill")]:
  model_path = os.path.join(models_path, y)
  with io.open(model_path, 'rb') as file:
      object_name = re.sub(r"\.", "_", y)
      object_name = re.sub(r"_dill", "", object_name)
      globals()[object_name] = dill.load(file)
      object_names.append(object_name)
#display(Markdown("Objetos cargados: \n\n>" + 
#         ", ".join(object_names)))

model_times = [o for o in object_names if o.endswith("_time")]
perf_table = pd.DataFrame({
  "Modelo": [re.sub("_time$", "", model_time) for model_time in model_times],
  "Tiempo": [globals()[model_time] for model_time in model_times]
})
df_n_params = pd.DataFrame(data={
    "Modelo": [x for x in model_n_params.keys()],
    "# Params": [x for x in model_n_params.values()]
})
perf_table = perf_table.merge(df_n_params, on="Modelo")
model_path = os.path.join(models_path, "model_n_params.dill")
model_n_params = dill.load(open(model_path, 'rb'))
model_histories = [o 
                   for o in object_names 
                   if o.endswith("_hist")
                  ]
model_metrics = [k 
                 for k in globals()[model_histories[0]].keys()
                 if re.search("^(val_|loss)", k) is None
                ]

for metric in model_metrics:
  perf_table["val_" + metric] = [np.mean(globals()[o]["val_" + metric]) for o in model_histories]
for metric in model_metrics:
  perf_table[metric] = [np.mean(globals()[o][metric]) for o in model_histories]
perf_table.rename({
  "mean_squared_logarithmic_error": "msle",
  "val_mean_squared_logarithmic_error": "val_msle"
}, axis=1, inplace=True)
perf_table.drop(["loss", "val_loss"], axis=1, inplace=True, errors='ignore')
perf_table.sort_values("val_mse", inplace=True)
perf_table["Tiempo"] = (perf_table["Tiempo"] // 60).astype('int').astype("str") + "m" + \
(perf_table["Tiempo"] % 60).round(3).apply(lambda x: f"{x:2.2f}") + "s"
perf_table.reset_index(inplace=True, drop=True)
#perf_table.round(4)
excluded_columns = ["iaqAccuracy", "datetime", "datetime-1", "delta", 
                    "imputated", "year"]
train, test = train_test_split(sinaica[[x 
                                        for x in sinaica.columns 
                                        if x not in excluded_columns]], 
                               train_size=0.8, random_state=175904, shuffle=False)
scaler_iaq = MinMaxScaler().fit(train[["IAQ"]])
perf_data_iaq = scaler_iaq.inverse_transform(perf_table.select_dtypes("float64"))
#scaler_gr = MinMaxScaler().fit(train[["gasResistance"]])
#perf_data_gr = scaler_gr.inverse_transform(perf_table.select_dtypes("float64"))
perf_data_iaq = pd.DataFrame(perf_data_iaq, 
                             columns=perf_table.select_dtypes("float64").columns)
#perf_data_gr  = pd.DataFrame(perf_data_gr, 
#                             columns=perf_table.select_dtypes("float64").columns)
perf_data_iaq.insert(0, "Tiempo", perf_table["Tiempo"], )
#perf_data_gr.insert(0,  "Tiempo", perf_table["Tiempo"], )
perf_data_iaq.insert(0, "Modelo", perf_table["Modelo"], )
#perf_data_gr.insert(0,  "Modelo", perf_table["Modelo"], )
perf_data_iaq["# Params"] = perf_table["# Params"]
#perf_data_gr["# Params"]  = perf_table["# Params"]
perf_data = perf_table.copy()
perf_data.sort_values("val_mse", inplace=True)
perf_data.reset_index(inplace=True, drop=True)
model_number_rows = [int(re.sub("[^0-9]", "", x)) for x in perf_table["Modelo"]]
#is_gr_row = [(x % 2) == 0 for x in model_number_rows]
is_iaq_row = [(x % 2) == 1 for x in model_number_rows]
#perf_data.iloc[is_gr_row] = perf_data_gr
perf_data.iloc[is_iaq_row] = perf_data_iaq
cols = ["Modelo", "Tiempo", "# Params", "val_mae", "mae"]
Markdown(perf_data.round(2)[cols].head(5).to_markdown())

	Modelo	Tiempo	# Params	val_mae	mae
0	model_best01b	9m34.02s	169,795	80.26	61.04
1	model_lstm01	2m24.19s	54,273	83.18	52.96
2	model_lstm03	6m13.69s	185,345	81.68	55.41
3	model_best01a	15m37.22s	575,745	83.38	53.35
4	model_conv01	2m53.16s	116,225	90.47	61.12

• Tabla Comparativa de los datos del sensor:

models = []
object_names = []
models_path = os.path.join(base_url, "models")

for y in [x for x in os.listdir(models_path) if x.endswith("dill")]:
  model_path = os.path.join(models_path, y)
  with io.open(model_path, 'rb') as file:
      object_name = re.sub(r"\.", "_", y)
      object_name = re.sub(r"_dill", "", object_name)
      globals()[object_name] = dill.load(file)
      object_names.append(object_name)
#display(Markdown("Objetos cargados: \n\n>" + 
#         ", ".join(object_names)))

model_times = [o for o in object_names if o.endswith("_time")]
perf_table = pd.DataFrame({
  "Modelo": [re.sub("_time$", "", model_time) for model_time in model_times],
  "Tiempo": [globals()[model_time] for model_time in model_times]
})
model_path = os.path.join(models_path, "model_n_params.dill")
model_n_params = dill.load(open(model_path, 'rb'))
#pd.DataFrame(model_n_params, columns=["# Parametros"])
df_n_params = pd.DataFrame(data={
    "Modelo": [x for x in model_n_params.keys()],
    "# Params": [x for x in model_n_params.values()]
})
perf_table = perf_table.merge(df_n_params, on="Modelo")
#df_n_params = perf_table.pop("# Params")
#perf_table.insert(2, "# Params", df_n_params)
model_histories = [o 
                   for o in object_names 
                   if o.endswith("_hist")
                  ]
model_metrics = [k 
                 for k in globals()[model_histories[0]].keys()
                 if re.search("^(val_|loss)", k) is None
                ]

for metric in model_metrics:
  perf_table["val_" + metric] = [np.mean(globals()[o]["val_" + metric]) for o in model_histories]
for metric in model_metrics:
  perf_table[metric] = [np.mean(globals()[o][metric]) for o in model_histories]
perf_table.rename({
  "mean_squared_logarithmic_error": "msle",
  "val_mean_squared_logarithmic_error": "val_msle"
}, axis=1, inplace=True)
perf_table.drop(["loss", "val_loss"], axis=1, inplace=True, errors='ignore')
perf_table.sort_values("val_mse", inplace=True)
perf_table["Tiempo"] = (perf_table["Tiempo"] // 60).astype('int').astype("str") + "m" + \
(perf_table["Tiempo"] % 60).round(3).apply(lambda x: f"{x:2.2f}") + "s"
perf_table.reset_index(inplace=True, drop=True)
#perf_table.round(4)
excluded_columns = ["iaqAccuracy", "datetime", "datetime-1", "delta", 
                    "imputated", "year"]
train, test = train_test_split(airdata[[x 
                                        for x in airdata.columns 
                                        if x not in excluded_columns]], 
                               train_size=0.8, random_state=175904, shuffle=False)
scaler_iaq = MinMaxScaler().fit(train[["IAQ"]])
perf_data_iaq = scaler_iaq.inverse_transform(perf_table.select_dtypes("float64"))
#scaler_gr = MinMaxScaler().fit(train[["gasResistance"]])
#perf_data_gr = scaler_gr.inverse_transform(perf_table.select_dtypes("float64"))
perf_data_iaq = pd.DataFrame(perf_data_iaq, 
                             columns=perf_table.select_dtypes("float64").columns)
#perf_data_gr  = pd.DataFrame(perf_data_gr, 
#                             columns=perf_table.select_dtypes("float64").columns)
perf_data_iaq.insert(0, "Tiempo", perf_table["Tiempo"], )
#perf_data_gr.insert(0,  "Tiempo", perf_table["Tiempo"], )
perf_data_iaq.insert(0, "Modelo", perf_table["Modelo"], )
perf_data_iaq["# Params"] = perf_table["# Params"]
#perf_data_gr["# Params"]  = perf_table["# Params"]
#perf_data_gr.insert(0,  "Modelo", perf_table["Modelo"], )
perf_data2 = perf_table.copy()
perf_data2.sort_values("val_mse", inplace=True)
perf_data2.reset_index(inplace=True, drop=True)
model_number_rows = [int(re.sub("[^0-9]", "", x)) for x in perf_table["Modelo"]]
#is_gr_row = [(x % 2) == 0 for x in model_number_rows]
is_iaq_row = [(x % 2) == 1 for x in model_number_rows]
#perf_data.iloc[is_gr_row] = perf_data_gr
perf_data2.iloc[is_iaq_row] = perf_data_iaq
cols = ["Modelo", "Tiempo", "# Params", "val_mae", "mae"]
perf_data2 = perf_data2.iloc[is_iaq_row]
Markdown(perf_data2.round(2)[cols].head(5).to_markdown())

	Modelo	Tiempo	# Params	val_mae	mae
1	model_dnn01	1m40.90s	4,609	74.06	61.15
2	model_best03a	14m0.58s	485,633	75.94	55.8
8	model_conv01	14m3.57s	294,401	127.78	5.79
9	model_conv03	6m33.58s	419,841	129.51	6.13
12	model_best03b	8m50.63s	162,115	164.43	9.28

Como se puede apreciar, incluso pudimos tener un mejor desempeño en los modelos que utilizan únicamente los datos del sensor y descartando los datos del gobierno. Lo cual resulta muy bueno porque no tendríamos que depender de los datos del gobierno. Además descartaríamos, de manera preliminar, que la contaminación exterior influye de manera significativa en la contaminación interior. Aunque destacamos que es necesario realizar más ensayos y pruebas de hipótesis para afirmarlo con una mayor certeza.

En la siguiente gráfica se puede notar que la red dense (DNN01) ofrece el mejor desempeño, sin embargo, destacamos que observamos resultados no estables; en algunas corridas de aprendizaje brindaba un excelente desempeño, y en otras no.

Modelos_list = perf_data2["Modelo"].tolist()
Modelos = pd.Categorical(perf_data2["Modelo"], 
                         categories=Modelos_list)
perf_data2["Modelo2"] = Modelos
perf_data["src"]  = "Sensor-Gov"
perf_data2["src"] = "Sensor"
_ = (
  ggplot(perf_data2.head(5), aes(x="Modelo2", y="val_mae", fill="Modelo2")) +
  geom_bar(stat="identity") +
  #geom_bar(aes(y="mae"), stat="identity") +
  labs(y="Error Medio\nAbsoluto", x="Modelos", 
       title="Gráfica Comparativa entre los Modelos\ncon Datos del Sensor"
      ) +
  theme(legend_position="none", axis_text_x=element_text(rotation=90), 
        axis_text_y=element_text(size=8),
        plot_title=element_text(size=11),) 
).draw()

Curvas de Desempeño

En la siguientes gráficas donde mostramos las curvas de desempeño con el error cuadrático medio (mse) en los 5 mejores modelos; intentando mostrar que el desempeño es mucho más estable en el modelo best01 que combina las redes convolucionales, recurrentes y densas comparado con el modelo denso simple.

models = []
object_names = []
for y in [x for x in os.listdir("models") if x.endswith("dill")]:
    with io.open(f"models/{y}", 'rb') as file:
        object_name = re.sub(r"\.", "_", y)
        object_name = re.sub(r"_dill", "", object_name)
        globals()[object_name] = dill.load(file)
        object_names.append(object_name)

models_path = "models"
model_path = os.path.join(models_path, "model_n_params.dill")
model_n_params = dill.load(open(model_path, 'rb'))

#Markdown("Objetos cargados: \n\n>" + 
#         ", ".join(object_names))
def models_plot(histories, a=None, b=None, 
                    metric='loss',
                    plot_validation=True):
  """
  Prints performance plot from a, to b on a list with object names in strings:
  i.e., ["model00", "model01", "modelltsm00", ..]
  
  Inputs:
  histories: an array with dicts containing the metrics key
  a: epoch start
  b. last epoch
  metric: plot this metric.
  plot_validation: boolean indicating if validation data should be plotted.
  a: from this epoch
  b: to this epoch    
  """
  
  """
  # Plot loss
  plt.subplot(imgrows, 2, 1)
  plt.title('Loss')
  plt.plot(history['loss'][a:b], label='Training', linewidth=2)
  if plot_validation:
    plt.plot(history['val_loss'][a:b], label='Validation', linewidth=2)
  plt.legend()
  plt.xlabel('Epoch')
  plt.ylabel(f'Loss')
  ### quantiles
  plt.xticks(ticks=quantiles-a,
              labels=quantiles)
  plt.grid(True)

  # Plot accuracy
  """
  models = [globals()[h] for h in histories]
  max_epochs = np.max([len(o[metric]) for o in models])
  if a is None:
      a = 0
  if b is None:
      b = max_epochs
  a = np.min((a,b))
  b = np.max((a,b))
  #print(f"a={a}, b={b}")
  
  imgrows = (len(histories)) / 2.
  imgrows = np.round(imgrows + .1, 0)
  imgrows = int(imgrows) #+ 1
  #print(f"length={len(histories)}")
  #print(f"imgrows={imgrows}")
  
  # x ticks
  quantiles = np.quantile(range(a, b), 
                          [.2, .4, .6, .8]).round(0).astype(int)
  quantiles = np.insert(quantiles, 0, [a])
  quantiles += 1
  quantiles = np.append(quantiles, [b-1])
  
  if plot_validation:
    model_order = {h: np.mean(v["val_" + metric]) for h, v in zip(histories, models)}
    #print(model_order)
  else:
    model_order = {h: np.mean(v[metric]) for h, v in zip(histories, models)}
  model_order = {m: model_order[m] for m in sorted(model_order, 
                                                   key=model_order.get)}
  
  # init the plot
  plt.figure(figsize=(14., 4.*imgrows), frameon=False, tight_layout=True)
  plt.suptitle("Gráfica de Desempeño de los Mejores Modelos con Datos del Sensor", 
               va='top', weight='bold', size='x-large')
  # create plot for every model
  #for i, (model_name, model) in enumerate(zip(histories, models)): 
  plt.autoscale(tight=True)
  for i, (model_name) in enumerate(model_order.keys()):
    model = globals()[model_name]
    
    plt.subplot(imgrows, 2, i+1)
    plt.title(model_name)
    if plot_validation:
      mean_label="Mean Validation"
    else:
      mean_label="Mean Training"
    plt.hlines(y=model_order[model_name], xmin=0, linestyles="dashed",
               xmax=b-a, label=mean_label, color="green")
    plt.plot(model[metric][a:b], label='Training', 
              linewidth=2)
    if plot_validation:
      plt.plot(model["val_" + metric][a:b], 
                label='Validation', linewidth=2)
      #print(f"{model_name}: {model_order[model_name]}")
    plt.legend()
    plt.xlabel('Epoch')
    plt.ylabel(metric)
    plt.xticks(ticks=quantiles-a, 
                labels=quantiles)
    plt.grid(True)

  plt.show()
  
model_histories = [o 
                   for o in object_names 
                   if o.endswith("_hist")
                  ]
model_metrics = [k 
                 for k in globals()[model_histories[0]].keys()
                 if re.search("^(val_|loss|mean_square_logarithmic_error)", k) is None
                ]
top_models = perf_data2.head(4)["Modelo"]
model_histories = [m + "_hist" for m in top_models]
models_plot(model_histories, metric='loss')

Conclusiones

Explicamos nuestras conclusiones por rubro:

• Análisis exploratorios sobre los datos:

  Encontramos sorprendente que el gobierno no se interese por publicar datos actualizados, por lo que fue necesario una buena inversión de tiempo para utilizarlos. Sin embargo, consideramos que los logramos superar con éxito mediante imputaciones y suplementar la información con otras estaciones de monitoreo ambiental como lo mencionamos en los puntos anteriores y en las siguientes líneas.

• Nuestras hipótesis y sus resultados:

  Descartamos nuestra hipótesis de la relación entre contaminación exterior e interior: no observamos una fuerterelación entre la contaminación exterior y la interior; pues tuvimos un mejor desempeño utilizando únicamente datos del sensor interior. Superando a los modelos con datos combinados del gobierno y el sensor. Esto nos resultó un poco sorprendente, pero encontramos y documentamos estos resultados a lo largo de este proyecto.

• Nuestras pruebas con modelos de redes neuronales profundas y el procesamiento de señales:

  Fueron necesarias muchísimas más pruebas de las que consideramos inicialmente, esto fue bueno, porque pusimos en práctica lo aprendido durante el curso, además de tener que buscar recursos en Internet. Esto fue enriquecedor y se retomará en los siguientes puntos.

• Los inconvenientes que encontramos y la influencia que ejercieron sobre el presente proyecto: describiremos a continuación las siguientes problemáticas a manera de resumen sobre lo mencionado en la sección Problemáticas.

  • Que no tuviéramos todos los datos disponibles en el sitio de la SINAICA nos complicó muchísimo no poder tener más datos.

  • El algoritmo propietario (cerrado) para calcular la variable IAQ a partir del sensor de la resistencia eléctrica en el aire hizo que perdiéramos tiempo en probar cuál de las variables era más viable predecir.

  • Aunque a caballo regalado no se le ve el diente: Los tiempos y capacidades de Google Colab no son los mejores, aún con Pro (que es bastante económico). Encontramos también que el desempeño computacional es aleatoria durante diferentes momentos.

  • Las fuentes externas no siempre son confiables, gratuitas o completas.

• Aprendizajes obtenidos durante el presente.

  • Nunca se debe subestimar la inversión de tiempo necesaria para limpiar, explorar, imputar, “corregir” y conocer los datos.

  • Cualquier mejora de desempeño conlleva un reto, a veces grande; tal como vimos en los miniproyectos y se comentó en clase.

  • No podemos confiar en que el gobierno o un ente externo entregue los datos a tiempo. Entendimos que es importante tener desde la fundación de un proyecto con al menos 70% de los datos disponibles, para empezar a trabajar con ellos.

  • Nunca subestimar las fuentes bibliográficas oficiales (en nuestro caso Keras y Tensorflow) como guías prácticas y bien escritas. Dado que muchas veces StackoOverflowow no tiene las mejores respuestas. [4], [5]

  • No tener miedo a aprender y a leer es vital para estar actualizado: pues muchas veces cambian las APIs, nuevos "approaches", nuevos paradigmas y nuevos descubrimiento de la inteligencia artificial.

  • Nos resultó muy útil tener un modelo baseline [5] que fue como una hipótesis nula, que nos permitió saber qué hacía un modelo vacío. Es bueno tener esta comparación. Consideramos esta idea muy valiosa.

  • Los detalles y el cuidado de estos son de suma importancia:

    i. Aprendimos a desescalar los datos para devolverlos a una escala real en las medidas originales.

    ii. Repetibilidad y reproducibilidad como principio científico: tuvimos cuidado en las semillas para ser lo más reproducible nuestro proceso, cuidamos que las posibles manipulaciones sobre el orden de los datos (como el train_test_split no hiciera shuffle ), escribimos siempre funciones que reutilizables a lo largo del código que permitieran repetir nuestros experimentos con los distintos modelos.

    iii. Es importante probar, probar y probar: alterando 1 solo detalle de la arquitectura, un hiperparámetro, una modificación de los datos, etc. Esto involucra hacer modelos constructivos, no hacerlos complejos desde un inicio.

• Siguientes pasos que nos vienen a la mente son:

  i. Tomar datos por al menos 1 año, para establecer más claramente la estacionalidad, incorporar (comprar) fuentes de datos del clima y contaminantes, tener más puntos de redundancia (no depender únicamente de 1 sensor en 1 sólo lugar).

  ii. Explorar mejores opciones para imputar datos secuenciales o series de tiempo.

  iii. Realizar más pruebas sobre las predicciones y ver cómo van las mediciones.

  iv. Hacer las pruebas con mayor capacidad de cómputo: tener ventanas más grandes en timeseries_dataset_from_array. Adicionalmente creemos pertinente probar la reescritura de esta función promediando, por ejemplo, a 1 min.

  v. Realizar más Hyperparameter Tuning (Learning Rate, combinar optimizadores en los modelos), buscar datos preentrenados, o entrenar la red con datos como el del ejemplo (del Departamento de Biogeoquímica del Instituto Max-Planck) [4] [5].

  vi. Hacer un tablero o que las predicciones sean predecibles por algún API.

  vii. Realizar predicciones con algún intervalo de confianza o grado de certidumbre de la predicción.

  viii. Mejorar el desempeño, que claramente no es el mejor.

  ix. Venderlo o publicarlo en algún journal.

Para concluir queremos también invitar al lector a que visite la página https://philwebsurfer.github.io/dlfinal/ donde hemos publicado todos los recursos necesarios para hacer reproducible y repetible este proyecto, así como los detalles de cada una de las secciones que fueron omitidos por espacio. Estos se han mencionado en las secciones anteriores de este artículo.

Bibliografía

• [1] Bosch Sensortec Community | How do I convert BME 680 gas resistance to IAQ?.

• [2] Bosch SensorTec Comunity | Solution to IAQ accuracy definition.

• [3] GitHub | Daniel Mancuso: Código fuente de OhmTech-io/uThingVOC Src/main.c

• [4] Keras contributors et al. Keras / Code examples / Timeseries / Timeseries forecasting for weather prediction. 2021.

• [5] Tensorflow Contributors. Tensorflow: Tutorial on Time series forecastingTime series forecasting. 2021.

• [6] Saadi S. M., et al. Development of indoor environmental index: Air quality index and thermal comfort index. 2017. doi:10.1063/1.4975276.

• Bosch BME680 Datasheet. 2021.

• Mancuso, Daniel. Indoor Air Quality Monitor | Hackster.io. 2019.

• Dirección de Monitoreo Atmosférico de la Secretaría del Medio Ambiente del Gobierno de la Ciudad de México.

• Sistema Nacional de Información de la Calidad del Aire del Gobierno Federal México.

• Pedregosa, Fabian; et al. Scikit-learn: Machine Learning in Python. 2021.

• Abadi, Martín, et al. TensorFlow: Large-scale machine learning on heterogeneous systems. 2015.

• McKinney, Wes. Data structures for statistical computing in python. 2010.

• Harris, Charles, et al. Array programming with NumPy. 2021. doi:10.1038/s41586-020-2649-2.

• Hunter, John. Matplotlib: A 2D Graphics Environment. 2007. doi:10.5281/zenodo.592536/zenodo.592536.

• Román-Rangel, Francisco. Notas y Código del Curso de Aprendizaje Profundo. 2021.

• González-Pérez, Felipe. Notas de aprendizaje de máquina. 2020.

• Mohd, Faizy. StackOverflow AI: How to use CNN for making predictions on non-image data?. 2021.

• Ondris, Ladislav. StackOverflow: How to apply CNN algorithm in python for non image dataset. 2020.

• Shady, Slim. Convolutional Neural Network With Tensorflow and Keras. 2021.

• Keras contributors et al. Keras API Reference: fit. 2021.

• Keras contributors et al. Keras API Reference: Convolution Layers > Convolution1D. 2021.

• Pedregosa, Fabian; et al. Scikit-learn: Machine Learning in Python. 2021.

• Abadi, Martín, et al. TensorFlow: Large-scale machine learning on heterogeneous systems. 2015.

• McKinney, Wes. Data structures for statistical computing in python. 2010.

• Harris, Charles, et al. Array programming with NumPy. 2021. doi:10.1038/s41586-020-2649-2

• Hunter, John. Matplotlib: A 2D Graphics Environment. 2007. doi:10.5281/zenodo.592536/zenodo.592536