library(tidyverse)
library(dplyr)
library(keras)Redes Neurais Recorrentes (RNN)
Uma Rede Neural Recorrente (RNR ou RNN, em inglês) é uma rede neural cujo objetivo é prever ou tirar cnclusões com base em entradas em sequências. A diferença entre estas redes e as Redes Neurais Artificiais (RNA ou ANN, em inglês)), é que as RNN consideram na sua arquitetura que a previsão da variável alvo no instante \(t\), \(y_t\), depende do seu valor em instantes anteriores (\(y_{t-k}\)). Já as ANN fazem previsões para a variável alvo \(y\) apenas considerando a observação de variáveis independentes.
O modelo matemático
Considere uma base de dados supervisionados. Vamos chamar de \(y_i\) a \(i\)-ésima observação da variável alvo. Vamos chamar de \(x_{i,j}\) a \(i\)-ésima observação da \(j\)-ésima variável independente. Em notação matricial, \(\mathbf{x}_i = (x_{i,1},x_{i,2},x_{i,2},\ldots,x_{i,m})\) representa o vetor de variáveis independentes da observação \(i\).
Vamos considerar que as observações ocorreram em instantes diferentes. Então, a observação \(i\) é referente ao instante \(i\). A frequência de observações pode ser qualquer, por exemplo, hora, dias ou meses.
Se usarmos uma RNA para prever o valor de \(y_i\) estamos suponto que \(y_i = f(\mathbf{x}_i)\), ou seja, que a variável alvo no instante \(i\) é função das variáveis independentes no instante \(i\). Nesse caso, uma possibilidade, considerando uma única camada oculta com 4 neurônios,
E o modelo matemático é:
\[ \begin{array}{lcl} y_i &=& g(w_{1,y}h_{i,1} + w_{2,y}h_{i,2} + w_{3,y}h_{i,3} + w_{4,y}h_{i,4} + \Theta_y)\\ h_{i,j} &=& g(w_{1,j}x_{i,1} + w_{2,j}x_{i,2} + w_{3,j}x_{i,3} + \Theta_j) \end{array} \]
Se usarmos uma RNR para prever o valor de \(y_i\) vamos supor que \(y_i = f(\mathbf{x}_{i},y_{i-1})\), ou seja, que a variável alvo no instante \(i\) é função não só das variáveis independentes no instante \(i\), como também é função de \(y_i\) em um instante anterior. Um exemplo semelhante ao anterior seria,
E o modelo matemático é:
\[ \begin{array}{lcl} y_i &=& g(w_{1,y}h_{i,1} + w_{2,y}h_{i,2} + w_{3,y}h_{i,3} + w_{4,y}h_{i,4} + \Theta_y)\\ h_{i,j} &=& g(w_{1,j}x_{i,1} + w_{2,j}x_{i,2} + w_{3,j}x_{i,3} + w_{h,j}h_{i-1,j} +\Theta_j) \end{array} \]
Mas para capturar ainda mais os padrões dos dados, podemos olhar para instantes maiores que o anterior modificando os dados de entrada da seguinte forma.
O modelo matemático também fica mais complexo.
\[ \begin{array}{lcl} y_i &=& g(w_{1,y}h_{i,1} + w_{2,y}h_{i,2} + w_{3,y}h_{i,3} + w_{4,y}h_{i,4} + \Theta_y)\\ h_{i,j} &=& g(w_{1,j}x_{i,1} + w_{2,j}x_{i,2} + w_{3,j}x_{i,3} + w_{h,j}h_{i-1,j} + w_{1,j}^{1}x_{i-1,1}+ w_{2,j}^{1}x_{i-1,2} + w_{3,j}^{1}x_{i-1,3} + w_{y,j}^{1}y_{i-1} + \Theta_j) \end{array} \]
Um exemplo prático
Considere a base Beijing PM2.5 Data, disponível em Link. Esta base contém medições horárias de poluição do ar em Pequim, especialmente da concentração de PM2.5, além de variáveis meteorológicas associadas.
| Nome da Variável | Papel | Tipo | Descrição | Unidade | Valores Ausentes |
|---|---|---|---|---|---|
No |
ID | Integer | Índice sequencial da observação | – | no |
year |
Feature | Integer | Ano da observação | year | no |
month |
Feature | Integer | Mês da observação | month | no |
day |
Feature | Integer | Dia do mês da observação | day | no |
hour |
Feature | Integer | Hora da observação (0–23) | hour | no |
pm2.5 |
Target | Integer | Concentração horária de partículas finas PM2.5 no ar | μg/m³ | yes |
DEWP |
Feature | Integer | Temperatura do ponto de orvalho (dew point) | °C | no |
TEMP |
Feature | Integer | Temperatura ambiente | °C | no |
PRES |
Feature | Integer | Pressão atmosférica | hPa | no |
cbwd |
Feature | Categorical | Direção combinada do vento (combined wind direction) | – | no |
Iws |
Feature | Continuous | Velocidade acumulada do vento (cumulated wind speed) | m/s | no |
Is |
Feature | Integer | Horas acumuladas de neve | hour | no |
Ir |
Feature | Integer | Horas acumuladas de chuva | hour | no |
Cada linha da base de dados passa a informação das variáveis acima em cada hora. Vamos considerar a variável alvo \(y_i\) como o valor de pm2.5 na linha \(i\), ou seja, a concentração de partículas finas PM2.5 no ar na faixa de hora definida pela linha \(i\). Já as variáveis independentes \(x_{i,j}\) são os valores de temperatura, pressão, … , tudo no instante \(i\).
Qual dos dois modelos parece mais apropriado?
\[ y_i = f(\mathbf{x}_i) \ \ \ \ \ \ \text{ou} \ \ \ \ \ \ y_i = f(\mathbf{x}_{i},\mathbf{x}_{i-1},\ldots,\mathbf{x}_{i-k}, y_{i-1},\ldots,y_{i-k}) \]
Aqui imaginamos que a a concentração de partículas finas PM2.5 no ar em uma certa hora depende não só das condições daquela hora como também das condições e da a concentração de partículas finas PM2.5 no ar nos instantes anteriores. Por isso vamos usar uma Rede Neural Recorrente.
Pacotes
Vamos primeiro carregar os pacotes necessários. O pacote keras será usado no treinamento da rede neural recorrente.
Leitura e organização da base de dados
base = read_csv("PRSA_data_2010.1.1-2014.12.31.csv")
glimpse(base)Rows: 43,824
Columns: 13
$ No <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 1…
$ year <dbl> 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010…
$ month <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
$ day <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
$ hour <dbl> 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18…
$ pm2.5 <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ DEWP <dbl> -21, -21, -21, -21, -20, -19, -19, -19, -19, -20, -19, -18, -19,…
$ TEMP <dbl> -11, -12, -11, -14, -12, -10, -9, -9, -9, -8, -7, -5, -5, -3, -2…
$ PRES <dbl> 1021, 1020, 1019, 1019, 1018, 1017, 1017, 1017, 1017, 1017, 1017…
$ cbwd <chr> "NW", "NW", "NW", "NW", "NW", "NW", "NW", "NW", "NW", "NW", "NW"…
$ Iws <dbl> 1.79, 4.92, 6.71, 9.84, 12.97, 16.10, 19.23, 21.02, 24.15, 27.28…
$ Is <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ Ir <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…base = na.omit(base)
base$time = as.POSIXct(
paste0(base$year,"/",base$month,"/",base$day," ",base$hour,":00"),
format = "%Y/%m/%d %H:%M"
)
base = base |> select(-c(year,month,day,hour))
glimpse(base)Rows: 41,757
Columns: 10
$ No <dbl> 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, …
$ pm2.5 <dbl> 129, 148, 159, 181, 138, 109, 105, 124, 120, 132, 140, 152, 148,…
$ DEWP <dbl> -16, -15, -11, -7, -7, -7, -7, -7, -8, -7, -7, -8, -8, -8, -9, -…
$ TEMP <dbl> -4, -4, -5, -5, -5, -6, -6, -5, -6, -5, -5, -5, -5, -5, -5, -5, …
$ PRES <dbl> 1020, 1020, 1021, 1022, 1022, 1022, 1023, 1024, 1024, 1025, 1026…
$ cbwd <chr> "SE", "SE", "SE", "SE", "SE", "SE", "SE", "SE", "SE", "SE", "SE"…
$ Iws <dbl> 1.79, 2.68, 3.57, 5.36, 6.25, 7.14, 8.93, 10.72, 12.51, 14.30, 1…
$ Is <dbl> 0, 0, 0, 1, 2, 3, 4, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1…
$ Ir <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ time <dttm> 2010-01-02 00:00:00, 2010-01-02 01:00:00, 2010-01-02 02:00:00, …summary(base) No pm2.5 DEWP TEMP
Min. : 25 Min. : 0.00 Min. :-40.00 Min. :-19.0
1st Qu.:11464 1st Qu.: 29.00 1st Qu.:-10.00 1st Qu.: 2.0
Median :22435 Median : 72.00 Median : 2.00 Median : 14.0
Mean :22279 Mean : 98.61 Mean : 1.75 Mean : 12.4
3rd Qu.:33262 3rd Qu.:137.00 3rd Qu.: 15.00 3rd Qu.: 23.0
Max. :43824 Max. :994.00 Max. : 28.00 Max. : 42.0
PRES cbwd Iws Is
Min. : 991 Length:41757 Min. : 0.45 Min. : 0.00000
1st Qu.:1008 Class :character 1st Qu.: 1.79 1st Qu.: 0.00000
Median :1016 Mode :character Median : 5.37 Median : 0.00000
Mean :1016 Mean : 23.87 Mean : 0.05534
3rd Qu.:1025 3rd Qu.: 21.91 3rd Qu.: 0.00000
Max. :1046 Max. :565.49 Max. :27.00000
Ir time
Min. : 0.0000 Min. :2010-01-02 00:00:00
1st Qu.: 0.0000 1st Qu.:2011-04-23 11:15:00
Median : 0.0000 Median :2012-07-23 17:30:00
Mean : 0.1949 Mean :2012-07-17 05:46:08
3rd Qu.: 0.0000 3rd Qu.:2013-10-17 20:15:00
Max. :36.0000 Max. :2014-12-31 23:00:00
NA's :5 Vejamos o comportamento de \(y\) nas 100 primeiras observações.
obs = 1:200
plot(x = obs,
y = base$pm2.5[obs],
type="b",las=2,xlab = "",ylab="cnt",axes=F)
axis(
1,
at = obs,
labels = base$time[obs],
las = 2,
cex.axis = 0.7
)
Definição da base de treino e teste
Quando temos dados em sequência a base de treino e teste não pode ser definida a partir de um sorteio aleatório. Precisamos estipular uma data de corte e os valores antes dela serão usados para treinamento e os valores depois para teste.
dim(base)[1] 41757 10t = 10000Vamos usar as 10000 primeiras observações para treino e as demais para teste.
Para facilitar, vou separar a base em três partes, variáveis independentes numéricas, categóticas e alvo.
base = base |> arrange(time)
X_n = base |> select(DEWP,TEMP,PRES,Iws,Is,Ir)
X_c = base |> select(cbwd)
Y = base$pm2.5base_treino = base[1:t,]
X_n_treino = X_n[1:t,]
X_c_treino = X_c[1:t,]
Y_treino = Y[1:t]Mudança de escala e tratamento das variáveis categóricas
Faremos a padronização já na base toda, mas lembrando que precisamos usar as medidas de média, desvio padrão, mínimo e máximo da base de treino.
Primeiro o tratamento nas variáveis independentes numéricas.
mediaX = colMeans(X_n_treino,na.rm = T)
sdX = sapply(X_n_treino, sd, na.rm = TRUE)
X_n = X_n |> apply(MARGIN = 1,FUN = function(l){(l - mediaX)/sdX}) |>
t() |> as.data.frame()Agora o tratamento na variável alvo.
minY = min(Y_treino)
maxY = max(Y_treino)
Y = (Y - minY)/(maxY - minY)Por fim, a transformação das variáveis categóricas em indicadoras.
X_c = model.matrix(~.,data = X_c)
X_c = X_c[,-1]
colnames(X_c)[1] "cbwdNE" "cbwdNW" "cbwdSE"Criar os dados de entrada da RNN
Vamos definir a variável janela, que fixa quantos passos para trás vamos usar na previsão dos passos à frente. O número de variáveis independentes será guardado em n_features.
janela <- 5
n_features <- ncol(X_n) + ncol(X_c)Os códigos a seguir constróem os dados como eles devem entrar na função de treinamento da rede neural recorrente. Vamos por partes. Primeiro será construída duas listas, X_seq e Y_seq, com os dados em sequência respeitando a janela estipulada.
X_seq <- list()
Y_seq <- list()
N = nrow(X_c)
for(i in 1:(N - janela)){
X_seq[[i]] <- cbind(
X_n[i:(i + janela - 1), ],
X_c[i:(i + janela - 1), ],
y=Y[i:(i + janela - 1)])
Y_seq[[i]] <- Y[(i+1):(i + janela)]
}X_seq[[1]] é uma matriz de dimensão janela \(times\) n_features que guarda os valores das variáveis independentes e também da variável alvo para as janela primeiras observações. Já Y_seq[[1]] é um vetor de tamanho janela que guarda os valores da variável alvo nos instantes seguintes.
X_seq[[1]] DEWP TEMP PRES Iws Is Ir cbwdNE
1 -0.9693659 -1.035932 0.1678180 -0.5021957 -0.08940524 -0.1316552 0
2 -0.9016553 -1.035932 0.1678180 -0.4873836 -0.08940524 -0.1316552 0
3 -0.6308127 -1.115083 0.2662132 -0.4725715 -0.08940524 -0.1316552 0
4 -0.3599701 -1.115083 0.3646084 -0.4427810 0.80197303 -0.1316552 0
5 -0.3599701 -1.115083 0.3646084 -0.4279689 1.69335130 -0.1316552 0
cbwdNW cbwdSE y
1 0 1 0.1307457
2 0 1 0.1501532
3 0 1 0.1613892
4 0 1 0.1838611
5 0 1 0.1399387Y_seq[[1]][1] 0.1501532 0.1613892 0.1838611 0.1399387 0.1103166A ideia é que vamos usar todas as informações em X_seq para prever o último valor de Y_seq.
Depois disso precisamos alocar cada um desses dois objetos em um array de dimensão 3. O array que vai guardar os valores de X_seq será chamado de X_array e tem dimensão (n, janela, n_features + 1), onde n é o número de elementos na lista X_seq.
n <- N-janela
X_array <- array(
NA_real_,
dim = c(n, janela, n_features + 1)
)
for(i in seq_len(n)) {
X_array[i, , ] <- as.matrix(X_seq[[i]])
}
dim(X_array)[1] 41752 5 10O array que vai guardar os valores de Y_seq será chamado de Y_array e tem dimensão (n, janela, 1), onde n é o número de elementos na lista Y_seq (igual ao número de elementos na lista X_seq).
Y_array <- array(
NA,
dim = c(n, janela,1)
)
for(i in 1:n){
Y_array[i, , 1] <- Y_seq[[i]]
}
dim(Y_array)[1] 41752 5 1A base de treino será construída com as primeiras 10000 matrizes desses objetos de 3 dimensões.
X_array_treino = X_array[1:t , ,]
Y_array_treino = Y_array[1:t , ,]Construíndo a rede neral em camadas
O pacote keras funciona um pouco diferente do que a gente fez até agora. Primeiro definimos a estrutura da rede neural. Cada objeto que entra na camada de entrada tem dimesão janela \(\times\) n_features + 1.
entrada <- layer_input(
shape = c(janela, n_features + 1)
)A saída será o resultado da entrada depois de passar pelas camadas ocultas. Aqui vamos supor duas camadas ocultas. A primeira com 5 neurônios e a segunda com 10.
saida <- entrada |>
layer_simple_rnn(
units = 5,
activation = "sigmoid",
return_sequences = TRUE
) |>
layer_simple_rnn(
units = 10,
activation = "sigmoid",
return_sequences = TRUE
) |>
layer_dense(
units = 1,
activation = "sigmoid"
)Em seguida a função keras_model define o modelo indicando quem serão as entradas e as saídas da rede neural.
modelo_rnn <- keras_model(
inputs = entrada,
outputs = saida
)
summary(modelo_rnn)Model: "model"
________________________________________________________________________________
Layer (type) Output Shape Param #
================================================================================
input_1 (InputLayer) [(None, 5, 10)] 0
simple_rnn_1 (SimpleRNN) (None, 5, 5) 80
simple_rnn (SimpleRNN) (None, 5, 10) 160
dense (Dense) (None, 5, 1) 11
================================================================================
Total params: 251 (1004.00 Byte)
Trainable params: 251 (1004.00 Byte)
Non-trainable params: 0 (0.00 Byte)
________________________________________________________________________________Ainda antes do treinamento, definimos como ele será realizado.
modelo_rnn |>
compile(
optimizer = optimizer_adam(
learning_rate = 0.001
),
loss = "mse",
metrics = "mae"
)Treinamento
O treinamento é então realizado com a função fit.
hist <- modelo_rnn |>
fit(
x = X_array_treino,
y = Y_array_treino,
epochs = 100,
validation_split = 0.2,
verbose = 0
)
plot(hist)
Previsão na base de treino
As previsões na base de treino.
pred <- modelo_rnn |> predict(X_array_treino)313/313 - 1s - 1s/epoch - 4ms/stepdim(pred)[1] 10000 5 1pred = pred[,janela,1]
pred = pred*(maxY - minY) + minYMedidas de qualidade na base de treino
Podemos usar como medida de qualidade na previsão o \(RMSE\)
\[ RMSE = \sqrt{\sum_{i=1}^n \dfrac{(y_i - \hat{y}_i)^2}{n}} \] ou um \(R^2\) adaptado. A adaptação será no setido de entender que no contexto de valores em sequência, o método mais simples e sem custo para calcular a previsão de \(y_i\) é considerar \(\hat{y}_{i-1}\).
\[ R^2 = 1 - \dfrac{\sum_{i=1}^n (y_i - \hat{y}_i)^2}{\sum_{i=1}^n (y_i - y_{i-1})^2} \]
Primeiro precisamos dos valores reais para comparação.
dim(Y_array_treino)[1] 10000 5Y = Y_array_treino[,janela]
Y = Y*(maxY - minY) + minYAgora já podemos calcular as medidas.
(rmse = sqrt(mean((Y - pred)^2)))[1] 26.04327sse = sum((Y - pred)^2)
Y1 = Y_array_treino[,janela-1]
sst = sum((Y - Y1)^2)
(R2 = 1 - sse/sst)[1] 0.9647332rmse = round(rmse,digits = 2)
R2 = round(R2,digits = 2)
fim = 100
plot(x = 1:fim, y=Y[1:fim],type="b",
main= "Previsão RNN - Treino - 1 passo a frente",
sub = paste0("RMSE = ",rmse," e R2 = ",R2),
xlab = "",ylab = "y")
points(x = 1:fim, y=pred[1:fim],type="b",col="red")
Previsão na base de teste
Definir primeiro a base de teste.
X_array_teste = X_array[(t+1):n , ,]
Y_array_teste = Y_array[(t+1):n , ,]pred <- modelo_rnn |> predict(X_array_teste)993/993 - 2s - 2s/epoch - 2ms/stepdim(pred)[1] 31752 5 1pred = pred[,janela,1]
pred = pred*(maxY-minY) + minYMedidas de qualidade na base de treino
A comparação com os valores reais.
dim(Y_array_teste)[1] 31752 5Y = Y_array_teste[,janela]
Y = Y*(maxY - minY) + minY(rmse = sqrt(mean((Y - pred)^2)))[1] 25.77611sse = sum((Y - pred)^2)
Y1 = Y_array_teste[,janela-1]
sst = sum((Y - Y1)^2)
(R2 = 1 - sse/sst)[1] 0.9627225rmse = round(rmse,digits = 2)
R2 = round(R2,digits = 2)
fim = 100
plot(x = 1:fim, y=Y[1:fim],type="b",
main= "Previsão RNN - Teste - 1 passo a frente",
sub = paste0("RMSE = ",rmse," e R2 = ",R2),
xlab = "",ylab = "y")
points(x = 1:fim, y=pred[1:fim],type="b",col="red")