Perspective Chapter: Airborne Pollution (PM2.5) Forecasting Using Long Short-Term Memory Deep Recurrent Neural Network Optimized by Gaussian Process

Abstract

Forecasting air pollution is a challenging problem today that requires special attention in large cities since they are home to millions of people who are at risk of respiratory diseases every day. At the same time, there has been exponential growth in the research and application of deep learning, which is useful to treat temporary data such as pollution levels, leaving aside the physical and chemical characteristics of the particles and only focusing on predicting the next levels of contamination. This work seeks to contribute to society by presenting a useful way to optimize recurrent neural networks of the short and long-term memory type through a statistical process (Gaussian processes) for the correct optimization of the processes.

1. Introduction

Recurrent neural networks (RNN), especially Long Short Term Memory (LSTM), have proved their efficiency in working on time-dependent values by (as its names indicate) the use of memory (sequences) gives enough information to the network to work properly finding patterns and trends in the values, which are not so obvious at first glance. Also, the gaussian processes are a useful statistical technique that allows the hyperparameters of the network since it has shown that the processing time is reduced and, at the same time, the accuracy may be improved [1]. Afterward, there is airborne pollution, which is a complex system that affects billions of people worldwide, especially in a metropolis such as Hotan, China. Shanghai, China, Ghaziabad, India, or in the case of this study, Mexico City, Mexico. Also, there we have a lot of types of particles interacting with each other in chemical, biological, and physical ways. The pollutants that are monitored by the SEDEMA’s network in the City of Mexico are nitrogen dioxide ($N{O}_2$), Ozone ($O_3$), sulfur dioxide ($S{O}_2$) and particulate matter ($P{M}_{2.5}$, and $P{M}_{10}$). Hence, we propose the use of an RNN-LSTM and optimizing its hyperparameters using Gaussian processes to increase the accuracy in the forecast of airborne pollution instead of the use of the current method.

2. Literature review

2.1 Airborne pollution

Air pollution has been a problem that has been increasing in recent decades mainly in large cities, bringing with it respiratory diseases [2, 3]. This contamination is accompanied by the same pollutant particles that have a useful life depending on physical parameters (size and shape) and their chemical composition. Mexico City has been studied for decades [4] due to its high levels of pollution that affect more than 20 million people. The sites used in this work are the following: Northeast (Gustavo A. Madero—GAM, FES Aragón—FAR, Xalostoc—XAL), Northwest (Tlalnepantla—TLA,), Center (Hospital General de México—HGM, Merced—MER), Southeast (Nezahualcóyotl—NEZ, Santiago Acahualtepec—SAC) and Southwest (Ajusco Medio—AJM, Pedregal—PED, Santa Fe, SFE). The map of the monitoring sites is shown in Figure 1.

Figure 1.

Location of the monitoring sites in Mexico City.

2.1.1 Multiple imputation by chained equations (MICE)

Dealing with the missing data problem often leads to two general approaches for imputing multivariate data: Joint modeling (JM) and fully conditional specification (FCS), also known as multivariate imputation by chained equations (MICE) [5]. It is known that is a JM-type problem when we must specify a multivariate distribution for the missing data and obtain the imputation of its conditional distributions through Markov Monte Carlo chains (MCMC) techniques. On the other hand, it is FCS, which specifies the multivariate imputation model on a variable-by-variable basis using a set of conditional densities, one for each incomplete variable. The imputation starts by iterating over the conditional densities, usually, a low number of iterations is enough. In order to explain the model, let us use the following notation [5]: Let $Y_j$ with ($j=1,\cdots$, p) be one of $p$ incomplete variables and $Y=$, $\cdots$, $Y_p$). The observed and missing parts of $Y_j$ are denoted by $Y_j^{\mathit{obs}}$ and $Y_j^{\mathit{mis}}$, respectively, then $Y^{\mathit{obs}}=\left(Y_1^{\mathit{obs}},\cdots ,Y_p^{\mathit{obs}}\right)$ and $Y^{\mathit{obs}}=\left(Y_1^{\mathit{mis}},\cdots ,Y_p^{\mathit{mis}}\right)$, these are the observed and missing data respectively in $Y$. The number of imputations is $m\ge 1$. The $h$-th imputed data set is denoted as $Y^{\left(h\right)}$ where $h=1,\cdots ,m$. Now let $Y_{-j}=\left(Y_1,\cdots ,Y_{j-1},Y_{j+1},\cdots ,Y_p\right)$ denote the collection of the $p-1$ variables in $Y$ except $Y_j$. Finally, let $Q$ denote the quantity of scientific interest. The mice algorithm has three main steps: imputation, analysis, and pooling. The analysis starts with an incomplete data set $Y_{\mathit{obs}}$. The second step is to compute $Q$ on each imputed data set, here the model is applied to $Y^{\left(1\right)},\cdots ,Y^{\left(m\right)}$ in the general identical. Finally, the third step is to pool the $m$ estimates $\widehat{Q^{\left(!\right)}},\cdots ,\widehat{Q^{\left(m\right)}}$ into one estimate $\overline{Q}$ and estimate its variance.

2.1.2 Recurrent neural networks

Recurring neural networks (better known as RNN) can be used for any type of data. In practical applications, the use of symbolic values is more common. In a recurrent neural network, there is a one-to-one correspondence between the layers of the network and specific positions in the sequence. The position in the sequence is also known as its timestamp. Finally, RNNs are complete Turing, which means that this type of network can simulate any algorithm with sufficient data and computational resources [6]. A representation of this kind of network is shown in Figure 2.

Figure 2.

Graphic representation of a recurrent neural network.

2.1.3 Long-short term memory (LSTM)

To represent the hidden states of the $\mathit{kth}$ hidden states (layer) the notation $h_t^{-\left(k\right)}$ is used and to simplify the notation it will be assumed that the input layer $\overline{x_t}$ can be denoted by $h_t^{-\left(0\right)}$ (this layer is not hidden) [7]. To obtain good results, a hidden vector of dimension p must also be included, which will be denoted by $\overline{c_t^{\left(k\right)}}$ and refers to the state of the cell. The state of the trap can be observed as the long-term memory within the network. The matrix that updates the values is denoted by $W^{\left(k\right)}$ and is used to permute the column vectors ${\left[h_t^{-\left(k-1\right)},h_{t-1}^{-\left(k\right)}\right]}^T$. The matrix that is obtained always results in dimensions $4p\times 2p$. A $2p$ size vector is then premultiplied by the $W^{\left(k\right)}$ matrix resulting in a $4p$ vector. Now to find the updates we have the following; for setting up intermediates Eq. 1, for selectively forget and add to long-term memory eq. 2, for selectively leak long-term memory to hidden state (Eq. 3).

$$\begin{array}{c}\mathit{\text{InputGate}}:\\ \mathit{\text{ForgetGate}}:\\ \begin{array}{c}\mathit{\text{OutputGate}}:\\ \mathit{\text{NewC}}-\mathit{\text{State}}:\end{array}\end{array}\left[\begin{array}{c}\overline{i}\\ \overline{f}\\ \begin{array}{c}\overline{o}\\ \overline{c}\end{array}\end{array}\right]=\left(\begin{array}{c}\mathit{\text{sigm}}\\ \begin{array}{c}\mathit{\text{sigm}}\\ \begin{array}{c}\mathit{\text{sigm}}\\ \mathit{\text{sigm}}\end{array}\end{array}\end{array}\right)W^{\left(k\right)}\left[\begin{array}{c}{\overline{h}}_t^{\left(k-1\right)}\\ {\overline{h}}_{t-1}^k\end{array}\right]$$E1

$$\overline{c_t^{\left(k\right)}}=\overline{f}\odot \overline{c_{t-1}^{\left(k\right)}}+\overline{i}\odot \overline{c}$$E2

$$\overline{h_t^{\left(k\right)}}=\overline{o}\odot \mathit{tanh}\left(\overline{c_t^{\left(k\right)}}\right)$$E3

Additionally, clarify that LSTM is an algorithm that belongs to recurrent neural networks or RNN [8]. The RNN’s refer to neural networks that take their previous state as input, this means that the neural network will have two inputs, the new information entered into the network and its previous state, which is shown in Figure 2. With this model we can have short-term memory in the neural network [9]. These neural networks have applications in sequential predictions, that is, predictions that depend on a temporal variable.

2.2 Bayesian optimization using Gaussian processes

2.2.1 Multidimensional Gaussian distribution

To talk about Gaussian processes, we must first define a multivariable Gaussian distribution in several dimensions. Formally, this distribution is expressed as in Eq. 4:

$$p\left(x\right)=\frac{1}{{\left(2\pi \right)}^{D/2}}{e}^{\frac{-1}{2}{\left(x-\mu \right)}^T}{\Sigma }^{-1}\left(x-\mu \right)$$E4

Where $D$ is the number of dimensions, $x$ is the variables, $\mu$ is the average vector, $\Sigma$ is the covariance matrix. Gaussian processes try to model a function f given a set of points [10]. Traditional nonlinear regression machine learning methods usually give a function that they think best fits these observations. But, there may be more than one function that fits the observations equally well. When we have more observation points, we use our posterior-anterior as our anterior, we use these new observations to update our posterior. This is the Gaussian process. A Gaussian process is a probability distribution over possible functions that fit a set of points. Because we have the probability distribution over all possible functions, we can calculate the means as the function and calculate the variance to show how confident we are when we make predictions using the function (Figure 3).

Figure 3.

Graphic example of a multidimensional gaussian distribution.

2.2.2 Gaussian process

Because we have the probability distribution over all possible functions, we can calculate the means as the function and calculate the variance to show how confident when predictions are made using the function as demonstrated by Wang [10], we must take into account that:

1. The (later) functions are updated with new observations.

2. The mean calculated by the posterior distribution of the possible functions is the function used for the regression.

The function is modeled by a multivariable Gaussian of the form shown in Eq. 5:

$$P\left(f|X\right)=N\left(f|\mu ,K\right)$$E5

Where $X=\left[x_{,}\cdots ,x_n\right]$, $f=$, $\cdots$, $f\left(x_n\right)$, $\mu =$, $\cdots ,$$m\left(x_n\right)$, $K_{i,j}=k\left(x_i,x_j\right)$. Being $X$ the points of the observed data, $m$ represents the average function and $K$ represents a definite positive kernel.

3. Materials and methods

3.1 Materials

The data used to train the model were obtained from the Atmospheric Monitoring System (SIMAT for its acronym in Spanish) database is conformed of four subsystems: RAMA, REDMA, REDMET, and REDDA, all are given by its website but we just focus on the Automatic Environmental Monitoring Network (RAMA for its acronym in Spanish).

3.2 Methodology

3.2.1 Data acquisition and preprocessing data

First, the database on air quality, RAMA (SEDEMA, 2021) [11] must be downloaded, which is public. This file (dataset) contains values captured by all stations capable of monitoring $P{M}_{2.5}$. In case the dataset contains more variables in addition to the one already mentioned, a preprocessing process must be carried out. First, the values of interest ($P{M}_{2.5}$. and air direction) should be classified, excluding any other. Once the data has been classified, it will be necessary to determine if there are missing data and in which cases it is convenient to impute because if the amount of data to be imputed exceeds 40% of the total data, it is advisable not to use that station since there is a loss very large data and a case of over-learning or data that does not reflect reality could be presented. To impute the missing data, the MICE algorithm will be used and once the algorithm has been applied, a new dataset will have to be generated with the imputed data (complete). Because the RNN-LSTM works by taking a tensor as input and already having an absent dataset of missing data, now it will be necessary to divide this data into three different datasets which would serve for training, testing, and data validation. To conclude with the preprocessing, the data will be normalized and later converted into tensors. To normalize the data, the min-max normalization will be used, which takes the minimum value of the data as “zero” and the maximum value as “one” and it is based on these that the normalization is performed. For tensors, they must take into account the batch value, which in turn takes values of $2^n$ with $\mathit{nϵN}$. The value of the sequence to use as a parameter should also be considered when creating the tensors.

3.2.2 Instantiate LSTM and optimize the model

To begin with the training and optimization of the network, we are going to start the network with values of the hyperparameters selected at random, this is only to make the network start and work since later the values of each selected hyperparameter will be rewritten in optimization until the optimal point is found within the search space that is established. By having the data imputed, divided, normalized, and transformed into tensors now, using the training dataset, as its name implies, we will begin the process of training the network, which will go hand in hand with the number of iterations (points) that have been assigned to the optimization process. In each iteration, an adjustment will be made in some of the hyperparameters, and saving the score of each model to end up with the best one.

3.2.3 Model evaluation

Having already a trained and optimized model, we can now determine the efficiency of the model by calculating the RMSE that the model has by comparing the predicted data against the real data. To generate predictions, we are going to use the evaluation set (Figure 4).

Figure 4.

Proposed methodology.

4. Results and discussion

There was a considerable amount of missing values all across the data. After applying the MICE algorithm the data full fill. In order to confirm that the imputation was successfully checked that the distribution of the data did not change as is shown in Figures 5 and 6. These processes will be repeated for each station. This process was repeated for all the stations which were selected getting similar results in all cases. Once the preprocessing was finished, the training was started and by optimizing the model, a search space was stabilized for each hyperparameter of the RNN-LSTM considered. Now, with the model ready and tested, the results of the optimization process can be seen within the network. As is shown in Figures 7 and 8, for the LSTM optimized by GP, the loss function during training decreases rapidly in each epoch until it converges and the change between epochs is no longer so noticeable compared with simple LSTM (Figures 9 and 10). In both cases when the converges are reached, it means that the network has already stopped learning. Finally, the validation set is used to estimate the skills of the network obtained in its training for the prediction of $P{M}_{2.5}$ levels. This is shown for the LSTM optimized by GP in Figures 11 and 12, and for the simple LSTM the results are shown in Figures 13 and 14. The results show are from AJM and PED stations.

Figure 5.

Data distribution of AJM station before and after being imputed.

Figure 6.

Data distribution of PED station before and after being imputed.

Figure 7.

Metrics obtained with LSTM model optimized by Gaussian process for the AJM station.

Figure 8.

Metrics obtained with LSTM model optimized by Gaussian process for the PED station.

Figure 9.

Metrics obtained with simple LSTM model for the AJM station.

Figure 10

Metrics obtained with simple LSTM model for the PED station.

Figure 11.

LSTM optimized by gaussian process, validation forecast for AJM station.

Figure 12.

LSTM optimized by gaussian process, validation forecast for PED station.

Figure 13.

Simple LSTM, validation forecast for PED station.

Figure 14.

Simple LSTM, validation forecast for PED station.

5. Conclusion

The potential of using Gaussian Processes to improve the hyperparameters tuning is based on a strong mathematical model compared with other methods of hyperparameter tuning; this gives more confidence when looking for an optimized deep learning model. On the other hand, to use this method requires more time to compute the Bayesian search and it may be a thing to consider.

About the project, it does not pretend to end the airborne pollution problem, but it introduces a useful tool for the government and the citizen in order to prevent and plane his day because the model can make predictions along more than three hundred hours with high accuracy, so this gives 12 days to the corresponding users to avoid or prevent the contamination levels.

Conflict of interest

The authors declare no conflict of interest

Consent for publication

Not applicable