Addressing model uncertainty in probabilistic forecasting using Monte Carlo dropout

Abstract

In recent years, deep learning models have been developed to address probabilistic forecasting tasks, assuming an implicit stochastic process that relates past observed values to uncertain future values. These models are capable of capturing the inherent uncertainty of the underlying process, but they ignore the model uncertainty that comes from the fact of not having infinite data. This work proposes addressing the model uncertainty problem using Monte Carlo dropout, a variational approach that assigns distributions to the weights of a neural network instead of simply using fixed values. This allows to easily adapt common deep learning models currently in use to produce better probabilistic forecasting estimates, in terms of their consideration of uncertainty. The proposal is validated for prediction intervals estimation on seven energy time series, using a popular probabilistic model called Mean Variance Estimation (MVE), as the deep model adapted using the technique.

Keywords

Probabilistic forecasting prediction interval monte carlo dropout mean variance estimation method epistemic uncertainty

1. Introduction

Several human activities depend on the forecasting of future values. A big ammount of research in machine learning has focused on how to perform these forecasting tasks using data of previous historical measurements. In general, these approaches are supervised, meaning that training data is used to build a model that transforms observations of values related to the forecasted phenomenon into predictions of the future values. These predictions, though, are uncertain, because of two reasons: on the one hand, it is impossible to know all the information related to the process generating the data. For example, when predicting load demand of a household, it depends on the activities of those who live in the house, and those activities depend on many factors that are not completly known, such as weather, their health, neighborhood activities, and more. This kind of uncertainty is called inherent or aleatoric uncertainty of the process, and it is present in most real life forecasting situations. On the other hand, many models can be used to produce predictions using the same previous information. Machine learning methods usually select one model from all possible ones by minimizing a desired loss function, chosen in advance, on training data. When doing so, it leaves out many other models that may fit the same data, and that, due to incomplete data, may be more close to the real process. This type of uncertainty, related to the model to use itself, is called epistemic uncertainty. In most probabilistic forecasting tasks, instead of a probabilistic distribution for future values, the problem is tackled producing prediction intervals, which give decision makers a good insight into how confident they can be of model predictions. This has been done for tasks that present a big inherent uncertainty, such as wind speed, wind power and electric consumer load forecasting, as mentioned in [1, 2]. Therefore, failing to consider all sources of uncertainty can lead to narrow prediction intervals, making systems that depend on these forecasts to be less reliable. An example of this can be found in wind power generation forecasting, where prediction intervals are used to decide which among the other power generators of the system have to work to satisfy the overall demand in case wind power is not enough. Consequently, narrow prediction intervals result in higher energy costs and a less reliable electrical grid [3].

During recent years, deep learning models composed of neural networks have attracted attention due to their capacity to model complex non linear relationships between past values and future predictions. To perform probabilistic forecasts, a popular technique is to apply the Mean Variance Estimation (MVE) method, which makes use of a neural network architecture with two outputs and is trained to predict the mean and variance of each future value [4]. This method allows the model to learn to predict the inherent uncertainty of the process and tackle prediction in heteroscedastic time series, that is, that have non-constant variance, and seems to work well in practice, as shown in [1, 2]. Unfortunately, the MVE method alone does not tackle the model uncertainty source, which is important for proabilistic forecasting. In order to tackle this uncertainty, other approaches have been explored recently that are based on variational inference, and have been applied on tasks that are different from probabilistic forecasting. These approaches replace network weights values by a posterior distribution learnt during training. One such approach is Monte Carlo (MC) dropout, a technique that randomly disables connections within the network, allowing it to learn to extract relevant features with redundancy during training, and then, during evaluation, can be seen as a simulation of an ensemble of networks that share parameters. To observe the effects that modelling both inherent and model uncertainty has compared to only modelling the former one, the results of a toy example are presented in Fig. 1. Here, we draw comparisons between a single MVE network, and an MVE network using MC dropout trained on synthetic one dimensional training data. Both models learn to predict a mean and a variance similar to the one used to generate the data in regions with training data, but they differ on how they make predictions in other regions. The MVE network cannot provide a reliable uncertainty estimation in other regions, because it does not consider the epistemic uncertainty, and this is specially noted in the uncertainty of the mean estimation. It can be seen that MC dropout estimates of the mean are more uncertain in regions without training data. For probabilistic forecasting tasks, this translates into wider and more accurate prediction intervals.

Figure 1.

Mean, variance and final prediction in black, for two models trained using the same synthetic training data in orange (grey when printed in black and white). First row: An MVE network captures changing variance in data, but does not consider model uncertainty in regions lacking training data. Second row: An MVE network using MC dropout captures both variance in data, and uncertainty of other regions, specially for the mean prediction.

The hypothesis of this work is that applying MC dropout to neural network models used for probabilistic forecasting tasks improves their performance due to incorporing the model uncertainty. In particular, this is experimentally supported for the estimation of prediction intervals using a Long Short-Term Memory (LSTM) network, comparing between a traditional MVE network and one with MC dropout for four wind speed, one wind power, one consumer electrical load and one utility power load forecasting tasks. This work is an extension of the research presented in [5]. The main differences are:

•

An additional wind power dataset and an additional consumer electrical load dataset have been added to the experimental evaluation.

•

A thorough review of the related work was done and included in an individual section.

•

There is a deeper quantitatively and qualitatively comparison between models.

•

There is a comparison between the proposed model and an ensemble of models, for one dataset.

The paper is structured as follows: Section 2 makes a review of related work, Section 3 introduces how prediction intervals are estimated, Section 4 details the proposal, Section 5describes the experiments to validate the proposal, and Section 6 makes final remarks and describe future research directions.

2. Related work

Initial probabilistic forecasting approaches based on neural networks, considered only homoscedastic time series [6], which is different to the approach of this work that can handle heteroscedastic time series. The works of [7, 8, 9] can handle heteroscedasticity using different optimization algorithms and objectives functions to build intervals. This work, instead, uses stochastic gradient optimization, which is simpler than those ad-hoc algorithms. Another work [10], also use a mix of complex models, which are harder to interpret and to optimize. Simpler approaches using similar optimization techniques than this work are [11, 12], which are based on Extreme Learning Machines (ELM), but do not use a recurrent model, as in this work. In particular, [12] considers the model uncertainty problem, as in this work, but using bootstrap, which is an ensemble technique, heavier than the proposed variational approach. At the same time, a complete solution using also a bootstrap ensemble of ELM networks to handle the model uncertainty problem is provided in [13], while this work handles it using MC dropout. The work of [14] improves these models using a recurrent network, as in this work, instead of a feed forward one or ELM, but it uses a complex optimization algorithm, and not stochastic gradient optimization as in this work. An interesting approach using deep architectures is [15], which considers an LSTM stacked to a CNN in order to produce point forecasts, and then the probabilistic part is modelled using clusters of similar data, although they do not consider the model uncertainty problem. Some works rely on input scenarios and output simulations to produce probabilistic forecasts, such as [16, 17, 18], using neural network models that are not essentially probabilistic, as in this work, and consider the model uncertainty using ensembles. In [19], a recurrent network is used as well, and its output is optimized to produce quantiles based on an adjustable distribution shaped by splines, which is flexible. In constrast, this work uses a simpler normal distribution, as the focus is to validate the model uncertainty issue. In the same line, [20] uses a single layer feed forward network and [21] uses a more complex neural network to produce quantiles, minimizing the pinball loss, while the output distribution is estimated using Kernel Density Estimation. Both works have networks with many outputs to produce quantiles, differently from this work where the network has only two ouputs: mean and variance. Another approach is [22], that uses a Deep Belief Network model, much harder to optimize and use than the simpler LSTM network used in this work. Closer to this work, [23] uses an MVE network, although not exploring the model uncertainty problem. Also, [24] uses a similar MVE network and addresses the model uncertainty problem, just as in this paper, but with two differences: it is trained using an ad-hoc metric for prediction intervals, which cannot be extended easily to other probabilistic forecasting tasks, and the model uncertainty problem is tackled with bootstrap, building an ensemble of models trained by resampling training data. Similarly, [25] addresses both sources of uncertainty using an MVE-like model, but the network architecture is more complex, having a wavelet transform and then applying a 2D CNN, which is different from the recurrent network used in this work. Finally, [26] is the closest approach to this work, using an MVE model for the inherent uncertainty and MC dropout to handle the model uncertainty, differing from this work because they use a GRU recurrent network instead of an LSTM, and they only apply it to a wind speed forecasting dataset, while in this work it is applied to seven datasets, including wind speed, wind power and load demand forecasting tasks.

A few works have tackled the model uncertainty problem without considering the inherent process uncertainty. For example, the ensemble model from [27] produces many point forecasts, changing the input for them, but each output stays deterministic, meaning the inherent process uncertainty is not considered, as in this work. Others assume homoscedasticity, that is, a fixed distribution for all times. In this scenario the use of ensemble is explored [28, 29, 30].

Aside from the forecasting community, there is a branch of research that explore neural network parameters uncertainty. In general, they make assumptions regarding the distribution of those parameters, and use variational inference to approximate their posterior distribution, as explained by [31]. One example is [32], which assumes Gaussian parameters distributions. Also, as explained by [33], many regularization techniques usual in the deep learning community can be interpreted as assumptions on the parameters distributions. In particular, MC dropout has been adapted to this framework and applied to many kinds of networks. This work is based on he application of MC dropout for recurrent neural networks, shown possible by [34]. Only few works have tackled the problem of model uncertainty using this variational approach for prediction intervals estimation. Among them, [35] deserves mention, which is applied for Uber data, focusing mostly on the coverage of produced intervals, and [5], which this work is an extended version of.

3. Prediction interval estimation

In this work, the probabilistic forecasting task being addressed is the estimation of prediction intervals. This task can be posed as the following problem: Given an univariate time-series $\{x_{t}\}$ , assumed to be generated by an underlying stochastic process, the goal is to estimate prediction intervals $I_{t}$ for time $t$ given by a lower bound $L_{t}$ and an upper bound $U_{t}$ , with a given nominal confidence level $c$ . This nominal confidence defines the expected probability of ocurrence of a value at time $t$ within the interval. In other words, the interval is built in order to satisfy that $p\left(x_{t}<L_{t}\right)\leqslant\frac{1-c}{2}$ , $p\left(x_{t}>U_{t}\right)\leqslant\frac{1-c}{2}$ and $p\left(L_{t}\leqslant x_{t}\leqslant U_{t}\right)\geqslant c$ . As the real distribution $p\left(x_{t}\right)$ is unknown, the goal of a probabilistic forecasting task is to model the distribution of future values based on previous observations of time series values, that is, the conditional distribution $p\left(x_{t+1}\mid x_{1},\ldots,x_{t}\right)$ . The number of past values can be very large, or not even all available, so a Markov assumption is made, meaning that the conditional distribution may be modelled using just the last $M$ values, finally working with a conditional model $p\left(x_{t+1}\mid\mathbf{x}_{t}\right)$ , where $\mathbf{x}_{t}=\left(x_{t-M+1},\ldots,x_{t}\right)$ . The relation between $x_{t+1}$ and $\mathbf{x}_{t}$ may be very complex, so it is assumed that it can be separated in two parts, as it is done by several works of probabilistic forecasting [4], [35]:

$\displaystyle x_{t+1}=f\left(\mathbf{x}_{t}\right)+n\left(\mathbf{x}_{t}\right),$ (1)

where $f\left(\mathbf{x}_{t}\right)$ is a deterministic part and $n\left(\mathbf{x}_{t}\right)$ is a random additive noise with mean 0. The assumption of a deterministic $f\left(\mathbf{x}_{t}\right)$ means that, in theory, there is a deterministic regressor $f$ that, given input values $\mathbf{x}_{t}$ , always have the same output, having no randomness, and leaving all the randomness to the noise. In this work, as it is done by most works, the noise $n\left(\mathbf{x}_{t}\right)$ is assumed normal, with a variance depending on $\mathbf{x}_{t}$ , allowing to rewrite the noise term as $n\left(\mathbf{x}_{t}\right)=\sigma\left(\mathbf{x}_{t}\right)\epsilon$ , where $\epsilon$ is a standard normal random variable, and $\sigma\left(\mathbf{x}_{t}\right)$ a scale factor function that is also deterministic, finally rewriting Eq. (1) as:

$\displaystyle x_{t+1}=f\left(\mathbf{x}_{t}\right)+\sigma\left(\mathbf{x}_{t}% \right)\epsilon,$ (2)

from where it can be seen that the distribution of $x_{t+1}$ , conditioned on $\mathbf{x}_{t}$ , is a normal distribution $\mathcal{N}\left(f\left(\mathbf{x}_{t}\right),\sigma\left(\mathbf{x}_{t}\right% )^{2}\right)$ . In the context of deep learning, the functions $f$ and $\sigma$ are the outputs of a neural network, which has a set of parameters, including the weights and biases. Calling $\theta$ this set of parameters, it can be seen that the functions $f$ and $\sigma$ are in fact parameterized by them, making the outputs of the network to depend on the inputs $\mathbf{x}_{t}$ and the parameters $\theta$ . Because of this, relation Eq. (2) also depends on $\theta$ , and can be written as:

$\displaystyle x_{t+1}=f\left(\mathbf{x}_{t},\theta\right)+\sigma\left(\mathbf{% x}_{t},\theta\right)\epsilon.$ (3)

Given training data $\mathcal{D}$ , most models are based on algorithms that find a good set of values for $\theta$ , assuming that training data is representative enough of the required learning task. This is not the case for time series like the ones in this work, as the time when data is collected may not provide enough samples to cover the entirety of the multi dimensional input space, specially if the input is comprised of many time steps $M$ . Because of this, the estimation of $\theta$ can be wrong. In order to see this more clearly, instead of a fixed value for $\theta$ , let’s assume that $\theta$ is a random variable. First, the conditional variance of $x_{t+1}$ is given by:

$\displaystyle\text{Var}\left[x_{t+1}\mid\mathbf{x}_{t}\right]=\text{Var}\left[% f\left(\mathbf{x}_{t},\theta\right)\right]+\text{Var}\left[\sigma\left(\mathbf% {x}_{t},\theta\right)\epsilon\right]=\text{Var}\left[f\left(\mathbf{x}_{t},% \theta\right)\right]+\mathbb{E}\left[\sigma\left(\mathbf{x}_{t},\theta\right)^% {2}\right].$ (4)

The first term is the variance of the mean estimation, as $\mathbb{E}\left[x_{t+1}\mid\mathbf{x}_{t}\right]=\mathbb{E}\left[f\left(% \mathbf{x}_{t},\theta\right)\right]$ , and the second is the expected variance given by the model. Equation (4) shows that, while the variance for a fixed $\theta$ is just $\sigma\left(\mathbf{x}_{t},\theta\right)^{2}$ , considering the model parameters as not completely certain may lead to a bigger variance of $x_{t+1}$ , and thus wider prediction intervals. Also, with this new assumption, the conditional distribution of $x_{t+1}$ is now given by

$\displaystyle p\left(x_{t+1}\mid\mathbf{x}_{t}\right)=\int_{\theta}p\left(x_{t% +1}\mid\mathbf{x}_{t},\theta\right)p\left(\theta\mid\mathcal{D}\right)d\theta% \text{,}$ (5)

where $p\left(\theta\mid\mathcal{D}\right)$ is called the posterior distribution of the parameters $\theta$ after training with data $\mathcal{D}$ , and $p\left(x_{t+1}\mid\mathbf{x}_{t},\theta\right)$ is given by the model when fixing $\theta$ . In this work, this last probability is assumed normal, having $x_{t+1}\mid\mathbf{x}_{t},\theta\sim\mathcal{N}\left(f\left(\mathbf{x}_{t},% \theta\right),\sigma\left(\mathbf{x}_{t},\theta\right)^{2}\right)$ .

Usually, neural network models based on Mean Variance Estimation (MVE) are trained to get a single value for its parameters, $\theta^{*}$ , using a training algorithm that maximizes the log-likelihood. Therefore, $\theta$ is replaced by $\theta^{*}$ in Eq. (3), reducing the integral of Eq. (5) to just $p\left(x_{t+1}\mid\mathbf{x}_{t}\right)=p\left(x_{t+1}\mid\mathbf{x}_{t},% \theta^{*}\right)$ , and the variance of Eq. (4) only to the second term: $\text{Var}\left[x_{t+1}\mid\mathbf{x}_{t}\right]=\sigma\left(\mathbf{x}_{t},% \theta^{*}\right)^{2}$ . It can be seen that the first term of the variance in Eq. (4) is completely ignored, making the corresponding prediction intervals to miss an important source of uncertainty. Ensemble models can be, and have been, used to tackle this problem, approximating the integral of Eq. (5) by a weighted sum of probabilities given by many models, each one with its own set of values for $\theta$ . While this approach offers a solution, in this work, another approach is preferred, which promises to be much simpler: Monte Carlo dropout, proposed by [36]. Dropout is a regularization technique, introduced by [37], that modifies the neural network each time it is evaluated, disabling some of its connections, or, in other words, making some weight values zero. The evaluation of the network can be seen as being over a subnetwork, composed by the active connections and units that remain after the deactivation. Evaluating the network many times over different subnetworks can be seen as equivalent to evaluating an ensemble of networks, thus sampling the posterior $p\left(\theta\mid\mathcal{D}\right)$ of Eq. (5) using another distribution implicitly defined by the selection of which connections are disabled. Hence, in this work, to create prediction intervals using this approach, the network parameters $\theta$ are sampled $K$ times, obtaining a set of different parameters $\{\theta_{1},\ldots,\theta_{K}\}$ . Then, the model in Eq. (3) is applied for each of them, obtaining $K$ possible future values, which are sorted and used to build the prediction interval, selecting the corresponding $\frac{1-c}{2}$ percentage value as its lower bound $L_{t+1}$ and the corresponding $\frac{1+c}{2}$ as its upper bound $U_{t+1}$ .

4. Proposal

The proposal of this work is to use a neural network model that combines the MVE method and MC dropout in order to build prediction intervals that include both inherent and model uncertainty, and that can handle heteroscedastic time series, allowing to produce better prediction intervals when presented with previously unseen data, compared to those obtained using the MVE method alone or including dropout only during training, as it is typically done.

Formally, the posterior distribution $p\left(\theta\mid\mathcal{D}\right)$ from Eq. (5) is unknown, hence it is approximated using another distribution $q\left(\theta\mid\phi\right)$ which is parameterized by a new set of parameters $\phi$ . The parameters $\phi$ are trained maximizing a lower bound of the train data log-likelihood, called evidence lower bound or variational bound, and given by:

$\displaystyle\text{ELBO}\left(\mathcal{D},\phi\right)=\mathbb{E}_{\theta\sim q% \left(\theta\mid\phi\right)}\left[\log p\left(\mathcal{D}\mid\theta\right)% \right]-D_{\text{KL}}\left(q\left(\theta\mid\phi\right)\mid\mid p\left(\theta% \right)\right)\leqslant\log p\left(\mathcal{D}\right).$ (6)

The inequality comes from Jensen’s inequality, and the fact that $\log$ is concave. The bound is composed by two terms: an approximate likelihood using distribution $q\left(\theta\mid\phi\right)$ instead of $p\left(\theta\mid\mathcal{D}\right)$ , and a divergence between an a priori assumption $p\left(\theta\right)$ and the posterior approximation. As explained by [36], using MC dropout is equivalent to assume implicit Bernoulli distributions for $q\left(\theta\mid\phi\right)$ , given by the probability of disabling a connection, while the second term is equivalent to $L_{2}$ , or weight decay regularization. Recalling from Eq. (5), a conditionally normal distribution was assumed for $x_{t+1}$ , given $\mathbf{x}_{t}$ and $\theta$ , so the log-likelihood term for each training data pair $\left(\mathbf{x}_{t},x_{t+1}\right)$ is given by

$\displaystyle\log p\left(x_{t+1}\mid\mathbf{x}_{t},\theta\right)=-\frac{1}{2}% \left[\log 2\pi+\log\sigma\left(\mathbf{x}_{t},\theta\right)^{2}+\frac{\left(x% _{t+1}-f\left(\mathbf{x}_{t},\theta\right)\right)^{2}}{\sigma\left(\mathbf{x}_% {t},\theta\right)^{2}}\right].$

Based on this likelihood, the variational bound has to be maximized, which can be treated as minimizing its negative version:

$\displaystyle\mathbb{E}_{\theta\sim q\left(\theta\right)}\left[\log\sigma\left% (\mathbf{x}_{t},\theta\right)^{2}+\frac{\left(x_{t+1}-f\left(\mathbf{x}_{t},% \theta\right)\right)^{2}}{\sigma\left(\mathbf{x}_{t},\theta\right)^{2}}\right]% +\text{Regularization}\left(\theta\right).$ (7)

This minimization is done using a stochastic gradient optimization algorithm, such as ADAM [38]. For each minibatch of size $B$ , the expectation is estimated by first sampling a value $\theta$ , say $\theta_{k}$ , which is then applied over the minibatch. In this way, Eq. (7) is estimated as:

$\displaystyle\sum_{b=1}^{B}\left[\log\sigma\left(\mathbf{x}_{t}^{\left(b\right% )},\theta_{k}\right)^{2}+\frac{\left(x_{t+1}^{\left(b\right)}-f\left(\mathbf{x% }_{t}^{\left(b\right)},\theta_{k}\right)\right)^{2}}{\sigma\left(\mathbf{x}_{t% }^{\left(b\right)},\theta_{k}\right)^{2}}\right]+\text{Regularization}\left(% \theta_{k}\right).$ (8)

To state clearly the steps performed by the training procedure, Algorithm 4 shows a step by step pseudocode version of it.

MVE Neural Network Model Training[1] Training time series $x_{1},\ldots x_{T}$ . Initialize network $h$ with randomly selected parameters $\theta$ . epoch $e=1$ $E$ mini-batch $k=1$ $K$ $\theta_{k}\leftarrow$ Sample weights applying MC dropout to current parameters $\theta$ . $L\leftarrow 0$ . Total mini-batch log-likelihood mini-batch sample $b=1$ $B$ $\textbf{x}\leftarrow\left(x_{t_{b}-M+1},\ldots,x_{t_{b}}\right)$ . Neural network input $f_{b},\sigma_{b}^{2}\leftarrow\text{Evaluate network }h\left(\textbf{x},\theta% _{k}\right)\text{to get mean and variance}$ . $L\leftarrow L-\frac{\log\sigma_{b}^{2}}{2}-\frac{\left(f_{b}-x_{t_{b}+1}\right% )^{2}}{2\sigma_{b}^{2}}$ . Compute log-likelihood $\text{Loss}\leftarrow-L+\text{Regularization}\left(\theta_{k}\right)$ . Update $\theta$ applying gradient descent algorithm, based on Eq. (8). Trained network parameters $\theta$ .

Note that, when no MC dropout is used, there is a similar estimation to be minimized by the optimization algorithm, although instead of using a sampled $\theta_{k}$ provided by dropout, a single current value for the parameters $\theta$ is used and optimized.

This strategy was implemented for a Long-Short Term Memory (LSTM), as explained in Fig. 2, which is a recurrent neural network that is commonly used for sequence processing, such as time series forecasting. To make a forecast, the network receives as input the $\mathbf{x}_{t}$ vector of $M$ previous values instead of the entire previous sequence, since adding more time steps in the past makes it slower and the LSTM is usually not able to learn to keep more information in its internal states despite making the input sequences longer. When applying dropout, the same connections are kept active or inactive for all time steps of the recurrence when evaluating one input vector $\mathbf{x}_{t}$ .

Table 1

Summary of the seven datasets used in this research to validate the proposal

Name	Domain	Resolution	Total	Mean	Min.	Max.	Missing (%)
B08	wind speed	1 hour	117215	6.72	0.00	17.0	0 (0%)
D05a	wind speed	1 hour	457155	7.20	0.00	26.1	0 (0%)
D08	wind speed	1 hour	320595	6.89	0.00	26.7	0 (0%)
E01	wind speed	1 hour	154266	7.78	0.00	25.5	0 (0%)
UCI individual household	electric power load	1 min.	2075259	1.09	0.08	11.1	25979 (1.25%)
GEFCom2014-L	electric power load	1 hour	61344	146.18	16.10	317.5	0 (0%)
Canela 1	wind power	1 hour	39889	3.13	0.00	19.0	1 (0.000025%)

Figure 2.

The input $\mathbf{x}_{t}$ is processed by the LSTM recurrence, producing outputs $f\left(\mathbf{x}_{t},\theta\right)$ and $\sigma\left(\mathbf{x}_{t},\theta\right)$ , which are intepreted as parameters of a normal distribution, which is then sampled from to produce a future value for $x_{t+1}$ . MC dropout modifies the network by deactivating some of the dashed connections, therefore using a sampled $\theta_{k}$ instead of a fixed $\theta$ . The same dropout mask is applied in all steps of the LSTM recurrence. Hence, different dropout masks will produce different $\tilde{x}_{t+1}$ samples, which are then used to build a predicion interval for $t+1$ .

4.1 Multi step forecasting

A model that considers model parameters uncertainty should have better estimates of uncertainty, specially in the long term, because errors usually accumulate through time steps. A model trained to perform one step ahead forecasts can be used to produce prediction intervals for horizons larger than a single time step in the following way: the output of a one step ahead prediction is injected as input for the model to produce the next one until the horizon is completed. This is done many times to simulate multiple future trajectories which can then be transformed into prediction intervals, as described in Fig. 3. This procedure is sometimes called rolling forward forecasting. Then, prediction intervals are built for each time step using the same sorting strategy explained for one step intervals.

Figure 3.

To build a trajectory of $\ell$ steps, a sample of $\theta$ is obtained using dropout, and then predictions are performed one step at a time, reinjecting the output as input for next time step, getting a trajectory of future values $\left(x_{t+1},\ldots,x_{t+\ell}\right)$ . After having $K$ trajectories, prediction intervals are built for each time step using gathered samples.

A step by step pseudocode procedure to produce prediction intervals for the following $\ell$ time steps is shown in Algorithm 3. It is based on getting a single value $\theta_{k}$ by selecting one dropout mask, then sampling a single value of $x_{t+1}$ from the normal distribution $\tilde{x}_{t+1}=f\left(\mathbf{x}_{t},\theta\right)+\sigma\left(\mathbf{x}_{t}% ,\theta\right)\epsilon$ , and then having this value re-injected as part of a new input $\textbf{x}_{t+1}=\left(x_{t-M+2},\ldots,\tilde{x}_{t+1}\right)$ to the model using the same dropout mask. A single iteration culminates when a trajectory of length $\ell$ is constructed, as described in Fig. 3. This process is repeated $K$ times with different dropout masks each, generating $K$ trajectories of sampled values, each with a different $\theta_{k}$ . At the end, $K$ sampled values will be available for each future time step, which are then used to build the prediction intervals corresponding to those time steps. For the MVE model without MC dropout, $K$ trajectories are computed as well, but removing the dropout step.

Multi-step Probabilistic Forecasting[1] Network with parameters $\theta$ , input values $x_{t-M+1},\ldots,x_{t}$ , number of future time steps $\ell$ . $k=1$ $K$ $\theta_{k}\leftarrow$ Sample weights applying MC dropout to $\theta$ $\textbf{x}\leftarrow\left(x_{t-M+1},\ldots,x_{t}\right)$ $j=1$ $\ell$ $f_{j}\leftarrow f\left(\textbf{x},\theta_{k}\right)$ First network output: mean $\sigma_{j}^{2}\leftarrow\sigma^{2}\left(\textbf{x},\theta_{k}\right)$ Second network output: variance $\epsilon\leftarrow$ Sample from $\mathcal{N}\left(0,1\right)$ $x_{t+j}^{k}\leftarrow f_{j}+\epsilon\sigma_{j}$ $\textbf{x}\leftarrow\left(x_{t-M+j+1}^{k},\ldots,x_{t+j}^{k}\right)$ with $x_{i}^{k}=x_{i}$ if $i\leqslant t$ $j=1$ $\ell$ $S\leftarrow$ Sort $\left[x_{t+j}^{1},\ldots,x_{t+j}^{K}\right]$ $L_{j}\leftarrow S\left[\frac{1-c}{2}K\right]$ $U_{j}\leftarrow S\left[\frac{1+c}{2}K\right]$ $\left[\left(L_{1},U_{1}\right),\left(L_{2},U_{2}\right),\ldots,\left(L_{\ell},% U_{\ell}\right)\right]$

5. Experiments

5.1 Data

To assess the proposal, seven datasets of four representative probabilistic forecasting tasks were considered, summarized in Table 1, showing different global distribution of values, as shown in Fig. 4. A brief description of each of them is given below:

•
B08, D05a, D08 and E01: These are four wind speed measurements datasets taken by prospecting towers located in different places of the Atacama Desert, in Northern Chile.1
¹
Data can be downloaded from http://walker.dgf.uchile.cl/Mediciones/.

Although the data collection is every 10 minutes, the datasets are provided in hourly averages. The distribution of values from these series varies, as it can be seen in Fig. 4, generally having a wide range of values.

Figure 4.
Histogram of time series values for all datasets.

•
UCI individual household electric power consumption: This dataset is composed of electrical load consumption measured every minute in a home for almost 4 years in France [39]. It presents many missing values, which were filled before the averaging process by copying the value at the same hour and minute from the previous day. For example, if the 10 am measurement of a day is missing, the value at 10 am of the previous day is used. The high resolution of this dataset makes this time series very hard to forecast, as the model would require information from many time steps in the past days for its input, in order to make a good minutely forecast. In this work, instead, the resolution was reduced to mean hourly values to be able to build models using the same architecture as the ones applied to the other datasets, taking the time series comprised by the hourly sum of each minute value as the input for the model. The distribution of values is highly two modal, as it can be seen in Fig. 4.
•
GEFCom2014-L: It corresponds to the power load of an electric utility in the US measured every hour for 7 years, and it was made available for a Global Energy Forecasting Competition [40]. Its distribution is smoother than the other series, with one dominant mode.
•
Canela 1: This dataset consists of hourly values of power injected to the electrical grid by a wind power plant located in Chile.2
²
Data can be downloaded from https://www.coordinador.cl/operacion/graficos/operacion-real/generacion-real-del-sistema/.

The distribution of this series is highly skewed, having many small values, and a more exponential distribution.

In order to evaluate the proposal over different regimes of a time series, each dataset was divided in five blocks, considering a small overlap between them, as shown in Fig. 5. For each block, the last 1000 measurements were left as test set and, from the remaining block data, 10% was considered as a validation set, used to select the best hyperparameters. For hyperparameters selection, a random search was done, as recommended by [41], and the selected parameters were those having the minimum average loss in the validation set, considering at least 5 repetitions with different random seeds. Hyper parameter values ranged according to values shown in Table 2.

Table 2
Range of hyper parameter values explored

Parameter name Min. Max. Scale

Hidden recurrent layers 1 3 Discrete

Hidden recurrent units 10 100 Discrete

Learning rate 0.00001 0.1 Log

Weight decay 0.0001 0.1 Log

Dropout rate 0.1 0.5 Linear

Epochs 25 45 Linear

The data was normalized for each block, considering only the training set to fix the minimum and maximum normalization parameters.

Figure 5.
Each time series is divided into five blocks of train, validation and test sets.

5.2 Metrics

Parameter name	Min.	Max.	Scale
Hidden recurrent layers	1	3	Discrete
Hidden recurrent units	10	100	Discrete
Learning rate	0.00001	0.1	Log
Weight decay	0.0001	0.1	Log
Dropout rate	0.1	0.5	Linear
Epochs	25	45	Linear

As the model proposed in this work was evaluated for the probabilistic forecasting task of prediction intervals estimation, the quality of the produced intervals has to be assessed considering the following three criteria:

•
Coverage: Each prediction interval is always associated to a Prediction Interval Nominal Confidence (PINC) that informs the associated probability of the next observation falling within the interval, as given by the model. Obviously, it is desired that future observations always fall inside the predicted intervals, but in reality this may not be the case. PINC can be seen as the expected percentage of future observations that may be inside of the predicted interval. In this work, all prediction intervals were evaluated for a nominal confidence (PINC) of 90%.
•
Sharpness: It would be easy to satisfy the coverage requirement by just making the intervals as wide as possible, but then they would be useless. Because of this, it is desired that the intervals be as narrow as possible, while still covering the required PINC.
•
Resolution: In this paper we work with heteroscedastic time series, which means that the width of the interval should change in time, adapting to the real uncertainty of each time period.

Many metrics exist in the literature that consider these features, differing on the importance given to each of them, and the penalization in case the coverage is not met by the intervals. In this work, the Winkler Loss, also called Interval Score (IS), was used as a general proxy of the quality of the interval. For each observation, it gives a base weight given by the width of the interval, and penalizes linearly the observations that fall out of the interval, according to:

$\displaystyle\text{Winkler}=\frac{1}{N}\sum_{i=1}^{N}\left[2\left(1-\text{PINC% }\right)\left(U_{i}-L_{i}\right)+4\begin{cases}L_{i}-y_{i}&\text{if }y_{i}<L_{% i}\\ 0&\text{if }L_{i}\leqslant y_{i}\leqslant U_{i}\\ y_{i}-U_{i}&\text{if }U_{i}<y_{i}\end{cases}\right],$

where $L_{i}$ and $U_{i}$ are the lower and upper bounds of the interval for time $i$ , and $y_{i}$ the test observation. Also, the following metrics were gathered, allowing for a detailed exploration of the characteristics of the intervals produced by the different models:

•
Prediction Interval Coverage Probability (PICP): It is a global measure of the percentage of observations falling inside the produced intervals: $\text{PICP}=\frac{1}{N}\sum_{i=1}^{N}\mathbb{I}\left(y_{i}\in\left[L_{i},U_{i}% \right]\right)$ . It is desired that PICP be always greater than PINC (90%).
•
Prediction Interval Normalized Average Width (PINAW): It measures the average width of produced intervals, as a percentage of the data range: $\text{PINAW}=\frac{1}{NR}\sum_{i=1}^{N}\left(U_{i}-L_{i}\right)$ , where $R=\max_{i\in\left\{1,\ldots,N\right\}}y_{i}-\min_{i\in\left\{1,\ldots,N\right% \}}y_{i}$ . As mentioned before, a PINAW as small as possible is desired, while having $\text{PICP}\geqslant\text{PINC}$ .
•
Coverage Width Criterion (CWC): It measures the quality of the interval globally, giving an exponential penalty when the coverage is not met. This metric is only low when the interval is working in terms of coverage, and then it favours narrower intervals: $\text{CWC}=\text{PINAW}\left(1+\mathbb{I}\left(\text{PICP}<\text{PINC}\right)e% ^{-\eta\left(\text{PICP}-\text{PINC}\right)}\right)$ . $\eta$ is a fixed value, usually selected as 50, which is also done in this work.

Figure 6.
Average Winkler, PICP and CWC metrics, considering different prediction horizons, for wind speed datasets. Standard deviation of each metric is displayed as shaded areas.

Figure 7.
Average Winkler, PICP and CWC metrics, considering different prediction horizons, for electrical load and wind power datasets. Standard deviation of each metric is displayed as shaded areas.

All MVE models, including the proposed model, may also be used to forecast in the traditional sense, that is, not probabilistically, producing the median or the most expected future value. For this reason, the models were also evaluated in terms of the following point prediction metrics:

•
Mean Absolute Error (MAE): It allows to compare the models in terms of their median point forecasting performance, $\text{MAE}=\frac{1}{N}\sum_{i=1}^{N}\left|y_{i}-\hat{y}_{i}\right|$ .
•
Root Mean Squared Error (RMSE): It allows to compare the models in terms of their mean point forecasting performance, $\text{RMSE}=\sqrt{\frac{1}{N}\sum_{i=1}^{N}\left(y_{i}-\hat{y}_{i}\right)^{2}}$ .

Table 3
Average metrics for each dataset and dropout mode (N: No dropout, Y: Regularization dropout, MC: Monte Carlo dropout), for 1 step ahead forecasting

Set Drop CWC PINAW (%) PICP (%) Winkler MAE RMSE

B08 N 19.517 (29.731) 30.2 (1.4) 86.1 (1.9) 8.07 (0.36) 1.60 (0.07) 2.02 (0.08)

Y 2.509 (1.787) 42.2 (2.2) 91.9 (1.2) 8.03 (0.37) 1.74 (0.06) 2.16 (0.07)

MC 0.959 (0.212) 35.4 (0.4) 91.1 (0.7) 7.14 (0.10) 1.58 (0.03) 1.98 (0.03)

D05a N 0.874 (0.425) 27.7 (1.2) 88.8 (1.3) 1.31 (0.02) 2.22 (0.11) 2.77 (0.13)

Y 0.414 (0.134) 27.7 (0.9) 90.6 (1.1) 1.22 (0.02) 2.13 (0.06) 2.66 (0.08)

MC 0.312 (0.038) 29.8 (0.8) 92.0 (0.7) 1.26 (0.03) 2.18 (0.10) 2.70 (0.10)

D08 N 1.219 (0.588) 25.5 (0.8) 88.5 (1.2) 1.29 (0.02) 2.26 (0.11) 2.82 (0.12)

Y 0.559 (0.307) 25.8 (1.3) 89.8 (1.4) 1.25 (0.03) 2.22 (0.04) 2.78 (0.05)

MC 0.344 (0.105) 30.3 (4.8) 92.1 (1.5) 1.34 (0.16) 2.28 (0.16) 2.83 (0.18)

E01 N 1.457 (0.442) 28.5 (0.6) 88.0 (0.7) 10.28 (0.22) 2.70 (0.12) 3.32 (0.17)

Y 0.356 (0.034) 34.3 (0.8) 92.1 (0.7) 9.62 (0.23) 2.74 (0.13) 3.34 (0.14)

MC 0.361 (0.035) 34.7 (0.5) 92.8 (0.7) 9.22 (0.19) 2.70 (0.11) 3.27 (0.11)

UCI N 0.968 (0.292) 36.6 (2.0) 90.7 (1.0) 28.68 (0.91) 33.7 (2.7) 43.6 (2.7)

Y 0.509 (0.182) 37.9 (1.3) 92.4 (1.0) 28.44 (0.67) 32.9 (2.0) 42.6 (1.7)

MC 0.424 (0.068) 40.5 (2.7) 93.4 (1.1) 28.65 (0.92) 35.7 (3.2) 45.0 (3.4)

GEF N 1.6 $\times$ 10 ${}^{9}$ (3.7 $\times$ 10 ${}^{9}$ ) 5.1 (0.5) 83.1 (7.4) 3.03 (0.72) 26.1 (10.2) 31.4 (11.4)

Y 0.081 (0.009) 8.1 (0.9) 95.2 (1.3) 3.37 (0.30) 20.7 (5.0) 26.1 (6.1)

MC 0.158 (0.047) 15.8 (4.7) 98.1 (0.6) 6.20 (1.97) 503.7 (1086) 11479.7 (25618)

Can. N 0.924 (0.411) 21.4 (1.4) 89.3 (1.3) 1.16 (0.03) 2.71 (0.56) 3.39 (0.56)

Y 0.261 (0.065) 23.4 (1.3) 92.4 (1.3) 1.12 (0.02) 2.36 (0.19) 3.00 (0.14)

MC 0.253 (0.009) 25.3 (0.9) 94.1 (0.7) 1.13 (0.02) 2.47 (0.11) 3.07 (0.10)

Table 4
Average metrics for each dataset and dropout mode (N: No dropout, Y: Regularization dropout, MC: Monte Carlo dropout), for 12 step ahead forecasting

Set Drop CWC PINAW (%) PICP (%) Winkler MAE RMSE

B08 N 28451 (47904) 42.3 (1.8) 82.5 (2.0) 13.32 (1.02) 1.61 (0.09) 2.04 (0.10)

Y 998 (656) 61.3 (3.3) 89.2 (2.2) 12.57 (0.79) 1.93 (0.08) 2.37 (0.09)

MC 598 (366) 45.4 (0.7) 86.3 (1.3) 11.74 (0.52) 1.58 (0.05) 1.99 (0.05)

D05a N 16.13 (23.83) 47.1 (2.3) 87.0 (2.0) 2.23 (0.11) 2.22 (0.11) 2.77 (0.13)

Y 2.19 (1.49) 46.8 (1.9) 89.0 (1.4) 2.11 (0.05) 2.15 (0.07) 2.68 (0.08)

MC 4.85 (6.13) 48.9 (2.1) 89.7 (1.5) 2.12 (0.08) 2.19 (0.09) 2.71 (0.10)

D08 N 348.5 (746.2) 44.8 (3.2) 85.4 (3.5) 2.35 (0.21) 2.31 (0.17) 2.89 (0.20)

Y 126.8 (186.5) 44.3 (2.6) 86.2 (2.4) 2.28 (0.14) 2.24 (0.06) 2.79 (0.07)

MC 12.6 (18.7) 47.7 (2.0) 88.7 (1.7) 2.19 (0.06) 2.25 (0.07) 2.79 (0.08)

E01 N 28.54 (31.42) 50.0 (1.8) 85.5 (2.0) 18.2 (1.6) 2.75 (0.15) 3.40 (0.17)

Y 18.72 (32.68) 53.1 (1.9) 87.5 (2.0) 16.6 (1.5) 2.72 (0.10) 3.31 (0.10)

MC 5.14 (7.16) 55.1 (1.3) 90.0 (1.6) 15.4 (1.1) 2.71 (0.09) 3.28 (0.11)

UCI N 2.49 (2.75) 47.5 (3.3) 90.9 (1.9) 35.9 (1.8) 34.3 (2.4) 44.0 (2.5)

Y 1.09 (0.67) 48.7 (3.4) 92.0 (1.3) 36.5 (1.8) 33.1 (2.2) 42.8 (1.9)

MC 1.08 (0.93) 51.7 (5.2) 92.7 (1.6) 36.7 (2.2) 35.8 (3.2) 45.0 (3.4)

GEF N 1.5 $\times$ 10 ${}^{17}$ (3.1 $\times$ 10 ${}^{17}$ ) 23.6 (2.0) 70.0 (12.8) 27.7 (16.8) 25.1 (8.8) 30.6 (9.9)

Y 5.3 $\times$ 10 ${}^{7}$ (1.2 $\times$ 10 ${}^{8}$ ) 31.3 (3.4) 84.2 (8.4) 19.1 (6.1) 21.3 (5.5) 26.7 (6.6)

MC 1085 (2421) 37.9 (3.9) 92.7 (3.8) 16.6 (2.0) 18.3 (4.0) 23.2 (5.3)

Can. N 8.6 $\times$ 10 ${}^{9}$ (1.9 $\times$ 10 ${}^{10}$ ) 44.8 (6.1) 85.4 (6.9) 2.66 (0.52) 2.62 (0.52) 3.31 (0.52)

Y 15732.43 (35176.26) 48.1 (6.5) 90.8 (3.4) 2.32 (0.30) 2.49 (0.35) 3.12 (0.35)

MC 1.02 (0.88) 51.5 (4.2) 93.2 (1.1) 2.23 (0.12) 2.50 (0.25) 3.10 (0.24)

5.3 Results

Set	Drop	CWC	PINAW (%)	PICP (%)	Winkler	MAE	RMSE
B08	N	19.517 (29.731)	30.2 (1.4)	86.1 (1.9)	8.07 (0.36)	1.60 (0.07)	2.02 (0.08)
	Y	2.509 (1.787)	42.2 (2.2)	91.9 (1.2)	8.03 (0.37)	1.74 (0.06)	2.16 (0.07)
	MC	0.959 (0.212)	35.4 (0.4)	91.1 (0.7)	7.14 (0.10)	1.58 (0.03)	1.98 (0.03)
D05a	N	0.874 (0.425)	27.7 (1.2)	88.8 (1.3)	1.31 (0.02)	2.22 (0.11)	2.77 (0.13)
	Y	0.414 (0.134)	27.7 (0.9)	90.6 (1.1)	1.22 (0.02)	2.13 (0.06)	2.66 (0.08)
	MC	0.312 (0.038)	29.8 (0.8)	92.0 (0.7)	1.26 (0.03)	2.18 (0.10)	2.70 (0.10)
D08	N	1.219 (0.588)	25.5 (0.8)	88.5 (1.2)	1.29 (0.02)	2.26 (0.11)	2.82 (0.12)
	Y	0.559 (0.307)	25.8 (1.3)	89.8 (1.4)	1.25 (0.03)	2.22 (0.04)	2.78 (0.05)
	MC	0.344 (0.105)	30.3 (4.8)	92.1 (1.5)	1.34 (0.16)	2.28 (0.16)	2.83 (0.18)
E01	N	1.457 (0.442)	28.5 (0.6)	88.0 (0.7)	10.28 (0.22)	2.70 (0.12)	3.32 (0.17)
	Y	0.356 (0.034)	34.3 (0.8)	92.1 (0.7)	9.62 (0.23)	2.74 (0.13)	3.34 (0.14)
	MC	0.361 (0.035)	34.7 (0.5)	92.8 (0.7)	9.22 (0.19)	2.70 (0.11)	3.27 (0.11)
UCI	N	0.968 (0.292)	36.6 (2.0)	90.7 (1.0)	28.68 (0.91)	33.7 (2.7)	43.6 (2.7)
	Y	0.509 (0.182)	37.9 (1.3)	92.4 (1.0)	28.44 (0.67)	32.9 (2.0)	42.6 (1.7)
	MC	0.424 (0.068)	40.5 (2.7)	93.4 (1.1)	28.65 (0.92)	35.7 (3.2)	45.0 (3.4)
GEF	N	1.6 $\times$ 10 ${}^{9}$ (3.7 $\times$ 10 ${}^{9}$ )	5.1 (0.5)	83.1 (7.4)	3.03 (0.72)	26.1 (10.2)	31.4 (11.4)
	Y	0.081 (0.009)	8.1 (0.9)	95.2 (1.3)	3.37 (0.30)	20.7 (5.0)	26.1 (6.1)
	MC	0.158 (0.047)	15.8 (4.7)	98.1 (0.6)	6.20 (1.97)	503.7 (1086)	11479.7 (25618)
Can.	N	0.924 (0.411)	21.4 (1.4)	89.3 (1.3)	1.16 (0.03)	2.71 (0.56)	3.39 (0.56)
	Y	0.261 (0.065)	23.4 (1.3)	92.4 (1.3)	1.12 (0.02)	2.36 (0.19)	3.00 (0.14)
	MC	0.253 (0.009)	25.3 (0.9)	94.1 (0.7)	1.13 (0.02)	2.47 (0.11)	3.07 (0.10)

Set	Drop	CWC	PINAW (%)	PICP (%)	Winkler	MAE	RMSE
B08	N	28451 (47904)	42.3 (1.8)	82.5 (2.0)	13.32 (1.02)	1.61 (0.09)	2.04 (0.10)
	Y	998 (656)	61.3 (3.3)	89.2 (2.2)	12.57 (0.79)	1.93 (0.08)	2.37 (0.09)
	MC	598 (366)	45.4 (0.7)	86.3 (1.3)	11.74 (0.52)	1.58 (0.05)	1.99 (0.05)
D05a	N	16.13 (23.83)	47.1 (2.3)	87.0 (2.0)	2.23 (0.11)	2.22 (0.11)	2.77 (0.13)
	Y	2.19 (1.49)	46.8 (1.9)	89.0 (1.4)	2.11 (0.05)	2.15 (0.07)	2.68 (0.08)
	MC	4.85 (6.13)	48.9 (2.1)	89.7 (1.5)	2.12 (0.08)	2.19 (0.09)	2.71 (0.10)
D08	N	348.5 (746.2)	44.8 (3.2)	85.4 (3.5)	2.35 (0.21)	2.31 (0.17)	2.89 (0.20)
	Y	126.8 (186.5)	44.3 (2.6)	86.2 (2.4)	2.28 (0.14)	2.24 (0.06)	2.79 (0.07)
	MC	12.6 (18.7)	47.7 (2.0)	88.7 (1.7)	2.19 (0.06)	2.25 (0.07)	2.79 (0.08)
E01	N	28.54 (31.42)	50.0 (1.8)	85.5 (2.0)	18.2 (1.6)	2.75 (0.15)	3.40 (0.17)
	Y	18.72 (32.68)	53.1 (1.9)	87.5 (2.0)	16.6 (1.5)	2.72 (0.10)	3.31 (0.10)
	MC	5.14 (7.16)	55.1 (1.3)	90.0 (1.6)	15.4 (1.1)	2.71 (0.09)	3.28 (0.11)
UCI	N	2.49 (2.75)	47.5 (3.3)	90.9 (1.9)	35.9 (1.8)	34.3 (2.4)	44.0 (2.5)
	Y	1.09 (0.67)	48.7 (3.4)	92.0 (1.3)	36.5 (1.8)	33.1 (2.2)	42.8 (1.9)
	MC	1.08 (0.93)	51.7 (5.2)	92.7 (1.6)	36.7 (2.2)	35.8 (3.2)	45.0 (3.4)
GEF	N	1.5 $\times$ 10 ${}^{17}$ (3.1 $\times$ 10 ${}^{17}$ )	23.6 (2.0)	70.0 (12.8)	27.7 (16.8)	25.1 (8.8)	30.6 (9.9)
	Y	5.3 $\times$ 10 ${}^{7}$ (1.2 $\times$ 10 ${}^{8}$ )	31.3 (3.4)	84.2 (8.4)	19.1 (6.1)	21.3 (5.5)	26.7 (6.6)
	MC	1085 (2421)	37.9 (3.9)	92.7 (3.8)	16.6 (2.0)	18.3 (4.0)	23.2 (5.3)
Can.	N	8.6 $\times$ 10 ${}^{9}$ (1.9 $\times$ 10 ${}^{10}$ )	44.8 (6.1)	85.4 (6.9)	2.66 (0.52)	2.62 (0.52)	3.31 (0.52)
	Y	15732.43 (35176.26)	48.1 (6.5)	90.8 (3.4)	2.32 (0.30)	2.49 (0.35)	3.12 (0.35)
	MC	1.02 (0.88)	51.5 (4.2)	93.2 (1.1)	2.23 (0.12)	2.50 (0.25)	3.10 (0.24)

In order to study the behavior of the proposal when using Monte Carlo dropout, it was compared against similar neural network models, first using regularization dropout only during training, and then without using dropout at all. While the models are trained to produce one step prediction intervals, they are run to produce multi step prediction intervals for 24 steps in the future, as explained previously, meaning that they produce intervals for each hour of the next day. Resulting average metrics for all time steps are displayed in Fig. 6, for wind speed datasets, and Fig. 7, for the other datasets. The best set of hyper parameters values was chosen for each dataset and time step shown.

It can be seen that, in general, MC dropout performs similarly to regularization dropout and no dropout for few time steps, but its handling of model uncertainty allows it to perform better in the long term, displaying reduced Winkler loss for many time steps. Looking at the PICP metric, it is clear that the coverage is higher using Monte Carlo dropout, translating into more observations falling inside the produced intervals, except for the D08 dataset. Interestingly, for this dataset, the Winkler loss is still low having a reduced PICP, meaning that, while having many observations out of the interval, its width was smaller than the linear penalization used by the metric. In general, as the coverage is still met for further time steps in the future, the CWC metric is also low in comparison with regularization or no dropout, even more than Winkler loss, because it penalizes exponentially when the coverage is not met. In the case of dataset B08, the reduced coverage also makes CWC higher and, for dataset D05a, CWC growths further for MC dropout, because unfortunately one of the repetitions had low coverage, which turned in an exponential penalty that, when averaged with the other repetitions, made it grow large. This effect is also seen on the big standard deviations related to CWC in general. For the same dataset D05a, and for the UCI dataset, the average Winkler loss was very similar for regularization and Monte Carlo dropout models, meaning that, while the coverage was increased, it was due to a bigger interval width. Interestingly, for datset UCI, Winkler and CWC losses were very similar for regularization and MC dropout. This could be related to the fact that the series has two modes instead of only one, making it difficul for models based on normal distributions to improve the produced intervals.

Table 5
Average metrics for each dataset and dropout mode (N: No dropout, Y: Regularization dropout, MC: Monte Carlo dropout), for 24 step ahead forecasting

Set	Drop	CWC	PINAW (%)	PICP (%)	Winkler	MAE	RMSE
B08	N	19131 (30056)	44.0 (1.9)	82.4 (2.0)	13.73 (1.08)	1.61 (0.08)	2.09 (0.10)
	Y	864 (555)	61.7 (3.2)	89.2 (2.1)	12.52 (0.78)	1.93 (0.08)	2.41 (0.09)
	MC	11718 (22479)	46.6 (1.0)	85.9 (1.4)	12.22 (0.54)	1.57 (0.04)	2.04 (0.04)
D05a	N	69.19 (86.98)	52.2 (2.7)	87.0 (2.0)	2.45 (0.12)	2.24 (0.09)	2.87 (0.11)
	Y	5.34 (5.76)	50.6 (2.8)	88.7 (1.9)	2.31 (0.09)	2.15 (0.07)	2.79 (0.09)
	MC	186.98 (405.49)	53.1 (2.6)	89.8 (1.8)	2.31 (0.12)	2.19 (0.09)	2.81 (0.11)
D08	N	96.5 (124.9)	49.5 (3.8)	85.7 (2.6)	2.53 (0.19)	2.33 (0.18)	2.99 (0.20)
	Y	258.5 (417.2)	47.3 (3.3)	85.9 (2.6)	2.43 (0.17)	2.24 (0.06)	2.88 (0.09)
	MC	10.3 (6.9)	50.7 (2.8)	88.4 (1.8)	2.36 (0.07)	2.27 (0.06)	2.91 (0.07)
E01	N	33.16 (47.38)	53.4 (2.5)	86.3 (2.2)	18.34 (1.71)	2.73 (0.15)	3.46 (0.19)
	Y	16.03 (30.55)	58.4 (3.0)	88.7 (2.1)	16.88 (1.42)	2.80 (0.12)	3.47 (0.13)
	MC	4.04 (5.68)	58.1 (1.2)	90.0 (1.6)	15.68 (0.97)	2.72 (0.09)	3.35 (0.10)
UCI	N	4.936 (7.938)	49.2 (5.8)	91.0 (2.2)	37.2 (2.7)	34.1 (3.3)	44.7 (3.0)
	Y	0.955 (0.500)	50.3 (3.7)	92.5 (1.3)	37.0 (1.7)	33.1 (2.2)	43.9 (1.9)
	MC	1.000 (0.814)	53.2 (5.2)	93.0 (1.6)	37.3 (2.3)	35.8 (3.2)	45.9 (3.4)
GEF	N	3.2 $\times$ 10 ${}^{17}$ (5.3 $\times$ 10 ${}^{17}$ )	30.4 (4.0)	67.9 (14.3)	39.4 (23.5)	24.3 (7.5)	30.8 (8.6)
	Y	5.0 $\times$ 10 ${}^{11}$ (9.9 $\times$ 10 ${}^{11}$ )	41.0 (6.5)	77.9 (14.4)	28.7 (12.4)	23.4 (7.2)	30.2 (8.7)
	MC	723510 (1.6 $\times$ 10 ${}^{6}$ )	44.1 (5.5)	90.8 (5.6)	19.8 (2.8)	18.0 (3.8)	23.6 (4.8)
Can.	N	2.1 $\times$ 10 ${}^{14}$ (4.7 $\times$ 10 ${}^{14}$ )	52.0 (11.5)	83.7 (8.9)	3.26 (1.08)	2.78 (0.72)	3.69 (0.68)
	Y	573329.47 (1.3 $\times$ 10 ${}^{6}$ )	54.9 (10.4)	90.0 (5.2)	2.68 (0.52)	2.52 (0.39)	3.36 (0.39)
	MC	1.27 (0.87)	58.5 (5.0)	93.3 (1.3)	2.54 (0.14)	2.51 (0.20)	3.27 (0.16)

Table 6

Average Winkler loss for GEFCom2014-L dataset, comparing dropout modes

Steps	No dropout	Regularization dropout	Monte carlo dropout
1	3.03 (0.72)	3.37 (0.30)	6.20 (1.97)
2	5.14 (1.67)	5.15 (0.50)	7.99 (1.84)
3	7.43 (2.92)	6.88 (0.79)	9.60 (1.85)
4	9.84 (4.40)	8.53 (1.17)	11.07 (2.05)
5	12.32 (6.02)	9.86 (1.56)	12.45 (2.37)
6	14.21 (7.02)	11.26 (2.02)	13.27 (2.09)
7	16.59 (8.62)	12.55 (2.50)	14.27 (2.21)
8	18.91 (10.23)	13.74 (2.97)	15.14 (2.34)
9	21.15 (11.84)	14.80 (3.40)	15.90 (2.46)
10	23.80 (13.94)	15.76 (3.78)	16.55 (2.57)
11	25.79 (15.41)	16.62 (4.14)	16.10 (1.85)
12	27.65 (16.75)	19.11 (6.07)	16.64 (2.01)
13	29.30 (17.62)	19.91 (6.45)	17.13 (2.17)
14	30.86 (18.68)	20.65 (6.81)	17.59 (2.32)
15	32.29 (19.60)	21.34 (7.16)	18.00 (2.46)
16	33.60 (20.42)	21.98 (7.51)	18.38 (2.60)
17	34.80 (21.17)	23.05 (9.25)	18.07 (2.14)
18	35.91 (21.88)	23.69 (9.75)	18.38 (2.24)
19	36.93 (22.56)	24.30 (10.24)	18.66 (2.34)
20	36.32 (21.12)	24.88 (10.71)	18.92 (2.44)
21	37.15 (21.71)	25.42 (11.16)	19.16 (2.53)
22	37.93 (22.31)	25.94 (11.60)	19.39 (2.62)
23	38.66 (22.90)	26.45 (12.03)	19.60 (2.71)
24	39.35 (23.48)	28.71 (12.43)	19.81 (2.80)

In order to compare the models in detail, Tables 3–5 show average metrics for 1, 12 and 24 steps ahead forecasting, including point forecasting metrics MAE and RMSE. While for 1 time step only dataset B08 presents a notable improvement, this is not the case for the other datasets, as seen in Table 3, meaning that considering the model uncertainty for one step ahead forecasting is not really necessary. Nevertheless, 12 step ahead results (Table 4) present mixed improvements, while 24 steps ahead results (Table 5) present a clear improvement using Monte Carlo dropout. This means that, although not necessary for short term prediction, considering model uncertainty using Monte Carlo dropout makes the model more robust for medium or long term forecasting tasks. This effect was specially notable for dataset GEFCom2014-L, an electrical load consumption dataset, which has the smoothest distribution, closer to the normal distribution used by the model. Also, as shown by point prediction metrics, adding Monte Carlo dropout does not reduce the performance of the model for point forecasting tasks, as shown by the MAE and RMSE metrics, meaning that models using Monte Carlo dropout can safely be used for multi step point forecasting tasks as well.

Figure 8.

Average metrics for Canela 1 dataset, considering different prediction horizons, and considering an ensemble of models. Standard deviation of each metric is displayed as shaded areas.

Table 7

Detailed metrics for selected time steps for Canela 1 dataset

Step	Drop	Winkler	PICP (%)	PINAW (%)	CWC	RMSE	Layers	Nodes
1	N	1.16 (0.03)	89.3 (1.3)	21.4 (1.4)	0.92 (0.41)	3.39 (0.56)	2.2	74.6
	Y	1.12 (0.02)	92.4 (1.3)	23.4 (1.3)	0.26 (0.06)	3.00 (0.14)	2.6	70.8
	MC	1.13 (0.02)	94.1 (0.7)	25.3 (0.9)	0.25 (0.01)	3.07 (0.10)	2.6	92.4
	E	1.10 (0.01)	94.2 (0.3)	25.1 (0.3)	0.25 (0.00)	3.13 (0.07)	2.0	51.2
3	N	1.60 (0.06)	88.3 (2.1)	30.1 (2.4)	58.79 (122.32)	3.38 (0.49)	2.0	61.8
	Y	1.52 (0.07)	92.1 (1.5)	31.7 (2.3)	0.52 (0.27)	3.20 (0.32)	2.2	75.4
	MC	1.50 (0.04)	93.6 (0.9)	33.6 (1.3)	0.36 (0.03)	2.98 (0.13)	2.2	96.2
	E	1.48 (0.01)	93.9 (0.3)	33.7 (0.5)	0.38 (0.03)	3.07 (0.05)	2.2	57.0
6	N	2.06 (0.23)	87.8 (3.7)	37.2 (3.6)	8.1 $\times$ 10 ${}^{5}$ (1.8 $\times$ 10 ${}^{6}$ )	3.26 (0.50)	2.2	77.2
	Y	1.89 (0.13)	91.5 (2.0)	39.8 (3.2)	8.67 (18.0)	3.23 (0.37)	2.0	59.6
	MC	1.87 (0.07)	93.3 (1.0)	41.7 (1.8)	0.48 (0.10)	3.03 (0.15)	2.0	94.8
	E	1.82 (0.01)	94.0 (0.4)	42.7 (0.8)	0.52 (0.05)	3.07 (0.05)	2.2	57.0
12	N	2.66 (0.52)	85.4 (6.9)	44.8 (6.1)	8.6 $\times$ 10 ${}^{9}$ (1.9 $\times$ 10 ${}^{10}$ )	3.31 (0.52)	1.8	72.6
	Y	2.32 (0.30)	90.8 (3.4)	48.1 (6.5)	1.6 $\times$ 10 ${}^{4}$ (3.5 $\times$ 10 ${}^{4}$ )	3.12 (0.35)	2.0	62.0
	MC	2.23 (0.12)	93.2 (1.1)	51.5 (4.2)	1.02 (0.88)	3.10 (0.24)	2.4	81.0
	E	2.21 (0.03)	93.9 (0.3)	52.6 (0.9)	0.67 (0.05)	3.06 (0.04)	2.0	55.2
24	N	3.26 (1.08)	83.7 (8.9)	52.0 (11.5)	2.1 $\times$ 10 ${}^{14}$ (4.7 $\times$ 10 ${}^{14}$ )	3.69 (0.68)	2.2	65.0
	Y	2.68 (0.52)	90.0 (5.2)	54.9 (10.4)	5.7 $\times$ 10 ${}^{5}$ (1.3 $\times$ 10 ${}^{6}$ )	3.36 (0.39)	2.6	68.8
	MC	2.54 (0.14)	93.3 (1.3)	58.5 (5.0)	1.27 (0.87)	3.27 (0.16)	2.6	70.4
	E	2.51 (0.05)	94.0 (0.3)	60.2 (2.1)	0.72 (0.06)	3.17 (0.05)	1.8	48.2

In Table 6, average Winkler loss for all time steps is displayed for dataset GEFCom2014-L. It is clear that using Monte Carlo dropout handles better the uncertainty for long term probabilistic forecasting tasks. Also, standard deviation, in brackets, grows more slowly when using Monte Carlo dropout than when not using it. This is in line with our consideration of model parameters uncertainty, because it means that the model is more robust to changes in the random seed that controls the randomized parts of the training, like initialization of parameters and mini-batch selection.

As mentioned before, applying MC dropout to a neural network may be interpreted as simulating an ensemble of subnetworks that share parameters. In order to compare the proposed model against an ensemble, similar experiments as described for each model were conducted using an ensemble of 20 similar neural networks using regularization dropout or no dropout at all, only for Canela 1 dataset, due to the very long time an ensemble takes to train. The results for all time horizons are displayed in Fig. 8, where it can be seen that the ensemble behaves better than the other models, specially for long term forecasting, as expected. The coverage (PICP) is increased for all time steps, slightly making the intervals wider (PINAW), while keeping metrics Winkler and CWC low for intervals, and low RMSE for point prediction. The standard deviation of all metrics was gathered considering 5 repetitions, and the ensemble model always had the lowest, as expected, because the ensemble considers all sources of uncertainty, while the MC dropout proposal is the one having the lowest standard deviation after the ensemble. This is inline with the assumption that using MC dropout considers model uncertainty, thus reducing the uncertainty of the training procedure and finally making the technique more robust. In Table 7 the same metrics are detailed for some selected time steps, and also the number of LSTM hidden layers and LSTM nodes of each layer. As MC dropout disables some connections, the selection of hyper parameters chose networks with more hidden nodes, thus making the network more similar to those of an ensemble of models after applying dropout.

6. Conclusion and future work

This work has explored how the behavior of a probabilistic forecasting recurrent deep learning model changes when considering the model uncertainty using a variational technique, in particular Monte Carlo dropout, for the task of producing prediction intervals. From the results, applying this technique for a recurrent model makes it to effectively handle the model uncertainty, and it is clear that this technique is specially important for multi-step prediction intervals generation, due to the increasing uncertainty of the future over the medium or long term, that is not completely handled by the base MVE neural network model.

From the results, it is concluded that, while the use of dropout is focused on handling the model uncertainty for probabilistic forecasting, its addition does not affect the point forecasting capabilities of the models, meaning that they can safely be improved for probabilistic forecasting using MC dropout, without interferring point forecasting metrics.

Additionally, the proposal seems to behave more similar to an ensemble of models for one dataset in terms of measured metrics. While this was expected, further comparison would be needed if their similarities and differences want to be clearly understood.

Although the proposal has been shown for univariate time series, the analysis for multivariate series is similar, and extending the use of MVE networks with MC dropout for them is straightforward, adding outputs of mean and variance estimates for each output dimension. Also, although in this work the deep architecture used is a recurrent one, the use of MC dropout can be done for other kind of networks, such as feedforward, convolutionals, extreme learning machines, etc.

As a future work, it is recommended to explore other variational techniques than dropout, in order to handle model parameters uncertainty, allowing more interpretability of the network parameters. Also, the assumption used in this work regarding the normal distribution of the outputs may be improved to handle time series with other distributions, specially those that are asymmetric or multimodal.

Footnotes

Acknowledgments

This work was supported in part by Conicyt doctoral scholarship 21170109, Fondecyt grant 1170123, Basal project FB0821, and UTFSM-PIIC grant 2019.

References

Kabir

H.M.D.

Khosravi

Hosen

M.A.

and Nahavandi

, Neural Network-Based Uncertainty Quantification: A Survey of Methodologies and Applications, IEEE Access 6 (2018), 36218–36234. doi: 10.1109/ACCESS.2018.2836917.

Khosravi

Nahavandi

Creighton

and Atiya

A.F.

, Lower upper bound estimation method for construction of neural network-based prediction intervals, IEEE Transactions on Neural Networks 22(3) (2011), 337–346. doi: 10.1109/TNN.2010.2096824.

Quan

Srinivasan

and Khosravi

, Incorporating Wind Power Forecast Uncertainties Into Stochastic Unit Commitment Using Neural Network-Based Prediction Intervals, IEEE Transactions on Neural Networks and Learning Systems 26(9) (2015), 2123–2135. doi: 10.1109/TNNLS.2014.2376696.

Nix

D.A.

and Weigend

A.S.

, Estimating the mean and variance of the target probability distribution, in: Proceedings of 1994 IEEE International Conference on Neural Networks (ICNN’94), Vol. 1, IEEE, Orlando, FL, 1994, pp. 55–60. ISBN 078031901X. doi: 10.1109/icnn.1994.374138.

Serpell

Araya

Valle

and Allende

, Probabilistic Forecasting Using Monte Carlo Dropout Neural Networks, in: Progress in Pattern Recognition, Image Analysis, Computer Vision, and Applications Nyström

Hernández Heredia

and Milián Núñez

, eds, Springer International Publishing, Havana, 2019, pp. 387–397.

Methaprayoon

Yingvivatanapong

Lee

W.-J.

and Liao

J.R.

, An Integration of ANN Wind Power Estimation Into Unit Commitment Considering, IEEE Transactions on Industry Applications 43(6) (2007), 1441–1448. http://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=4385004.

Khosravi

Nahavandi

Creighton

and Atiya

A.F.

, Comprehensive review of neural network-based prediction intervals and new advances, IEEE Transactions on Neural Networks 22(9) (2011), 1341–1356. doi: 10.1109/TNN.2011.2162110.

Quan

Srinivasan

and Khosravi

, Short-term load and wind power forecasting using neural network-based prediction intervals, IEEE Transactions on Neural Networks and Learning Systems 25(2) (2014), 303–315. doi: 10.1109/TNNLS.2013.2276053.

Shi

Liang

and Dinavahi

, Wavelet neural network based multiobjective interval prediction for short-term wind speed, IEEE Access 6 (2018), 63352–63365. doi: 10.1109/ACCESS.2018.2877735.

10.

Haque

A.U.

Nehrir

M.H.

and Mandal

, A hybrid intelligent model for deterministic and quantile regression approach for probabilistic wind power forecasting, IEEE Transactions on Power Systems 29(4) (2014), 1663–1672. doi: 10.1109/TPWRS.2014.2299801.

11.

Wang

Chen

and Chen

, Improved Chicken Swarm Algorithms Based on Chaos Theory and Its Application in Wind Power Interval Prediction, Mathematical Problems in Engineering 2019 (2019). doi: 10.1155/2019/1240717.

12.

Mahmoud

Dong

Z.Y.

and Ma

, An advanced approach for optimal wind power generation prediction intervals by using self-adaptive evolutionary extreme learning machine, Renewable Energy 126 (2018), 254–269. doi: 10.1016/j.renene.2018.03.035.

13.

Wan

Pinson

Dong

Z.Y.

and Wong

K.P.

, Probabilistic forecasting of wind power generation using extreme learning machine, IEEE Transactions on Power Systems 29(3) (2014), 1033–1044. ISBN 0885-8950. doi: 10.1109/TPWRS.2013.2287871.

14.

Shi

Liang

and Dinavahi

, Direct Interval Forecast of Uncertain Wind Power Based on Recurrent Neural Networks, IEEE Transactions on Sustainable Energy 9(3) (2018), 1177–1187. doi: 10.1109/TSTE.2017.2774195.

15.

Chen

Qiao

and Lu

, Probabilistic short-term wind power forecasting based on deep neural networks, in: 2016 International Conference on Probabilistic Methods Applied to Power Systems (PMAPS), IEEE, Beijing, 2016, pp. 1–8. ISBN 9781509019700. doi: 10.1109/PMAPS.2016.7764155.

16.

Gan

Wang

Yang

and Kang

, Embedding based quantile regression neural network for probabilistic load forecasting, Journal of Modern Power Systems and Clean Energy 6(2) (2018), 244–254. doi: 10.1007/s40565-018-0380-x.

17.

Xie

and Hong

, GEFCom2014 probabilistic electric load forecasting: An integrated solution with forecast combination and residual simulation, International Journal of Forecasting 32(3) (2016), 1012–1016. doi: 10.1016/j.ijforecast..

18.

Xie

and Hong

, Temperature Scenario Generation for Probabilistic Load Forecasting, IEEE Transactions on Smart Grid 9(3) (2018), 1680–1687.

19.

Gasthaus

Benidis

Wang

Rangapuram

S.S.

Salinas

Flunkert

and Januschowski

, Probabilistic Forecasting with Spline Quantile Function RNNs, Proceedings of Machine Learning Research 89 (2019), 1901–1910. http://proceedings.mlr.press/v89/gasthaus19a.html.

20.

and Li

, Probability density forecasting of wind power using quantile regression neural network and kernel density estimation, Energy Conversion and Management 164(January) (2018), 374–384. doi: 10.1016/j.enconman.2018.03.010.

21.

Hatalis

Lamadrid

A.J.

Scheinberg

and Kishore

, Smooth Pinball Neural Network for Probabilistic Forecasting of Wind Power, 2017, pp. 1–8. http://arxiv.org/abs/1710.01720.

22.

Wang

H.-Z.

Wang

G.-B.

G.-Q.

Peng

J.-C.

and Liu

Y.-T.

, Deep belief network based deterministic and probabilistic wind speed forecasting approach, Applied Energy 182 (2016), 80–93. doi: 10.1016/j.apenergy.2016.08.108.

23.

Khosravi

and Nahavandi

, An optimized mean variance estimation method for uncertainty quantification of wind power forecasts, International Journal of Electrical Power and Energy Systems 61 (2014), 446–454. doi: 10.1016/j.ijepes.2014.03.060.

24.

Khosravi

Nahavandi

and Creighton

, Prediction intervals for short-term wind farm power generation forecasts, IEEE Transactions on Sustainable Energy 4(3) (2013), 602–610. doi: 10.1109/TSTE.2012.2232944.

25.

Wang

H.-Z.

G.-Q.

Wang

G.-B.

Peng

J.-C.

Jiang

and Liu

Y.-T.

, Deep learning based ensemble approach for probabilistic wind power forecasting, Applied Energy 188 (2017), 56–70. ISBN 0306-2619. doi: 10.1016/j.apenergy.2016.11.111.

26.

Shahroudi

, Probabilistic Forecasting with Monte-Carlo Dropout in Neural Networks, Master’s thesis, University of Tartu, 2019.

27.

Liu

Nowotarski

Hong

and Weron

, Probabilistic load forecasting via quantile regression averaging of independent expert forecasts, IEEE Transactions on Smart Grid 8(2) (2017), 730–737. ISBN 2013112203. doi: 10.1109/TSG.2015.2437877.

28.

Sideratos

and Hatziargyriou

N.D.

, Probabilistic wind power forecasting using radial basis function neural networks, IEEE Transactions on Power Systems 27(4) (2012), 1788–1796. doi: 10.1109/TPWRS.2012.2187803.

29.

Cheng

Zang

Ding

Sun

Wang

Wei

and Sun

, Ensemble recurrent neural network based probabilistic wind speed forecasting approach, Energies 11(8) (2018). doi: 10.3390/en11081958.

30.

Wang

Zhang

Chen

Kirschen

D.S.

and Xia

, Data-Driven Probabilistic Net Load Forecasting With High Penetration of Behind-the-Meter PV, IEEE Transactions on Power Systems 33(3) (2018), 3255–3264. doi: 10.1109/TPWRS.2017.2762599.

31.

Graves

, Practical Variational Inference for Neural Networks, in: Advances in Neural Information Processing Systems 24 (NIPS 2011) Shawe-Taylor

Zemel

R.S.

Bartlett

P.L.

Pereira

and Weinberger

K.Q.

, eds, Curran Associates Inc., Granada, 2011, pp. 2348–2356. ISBN 9781618395993. https://papers.nips.cc/paper/4329-practical-variational-inference-for-neural-networks.pdf.

32.

Blundell

Cornebise

Kavukcuoglu

and Wierstra

, Weight Uncertainty in Neural Networks, in: Proceedings of the 32nd International Conference on Machine Learning (ICML 2015), Vol. 37 Bach

and Blei

, eds, PMLR, Lille, 2015, pp. 1613–1622. ISSN 1938-7228. ISBN 9781510810587. http://arxiv.org/abs/1711.01297.

33.

Gal

and Ghahramani

, On Modern Deep Learning and Variational Inference, in: Advances in Approximate Bayesian Inference workshop, Neural Information Processing Systems (NIPS), Vol. 2, 2015, pp. 1–9.

34.

Gal

and Ghahramani

, Bayesian Convolutional Neural Networks with Bernoulli Approximate Variational Inference, in: 4th International Conference on Learning Representations (ICLR 2016), San Juan, 2016, pp. 1–12.

35.

Zhu

and Laptev

, Deep and Confident Prediction for Time Series at Uber, in: 2017 IEEE International Conference on Data Mining Workshops (ICDMW), New Orleans, LA, 2017, pp. 103–110. ISSN 23759259. ISBN 9781538614808. doi: 10.1109/ICDMW.2017.19.

36.

Gal

and Ghahramani

, A Theoretically Grounded Application of Dropout in Recurrent Neural Networks, in: Advances in Neural Information Processing Systems 29 (NIPS 2016) Lee

D.D.

Sugiyama

Luxburg

U.V.

Guyon

and Garnett

, eds, Curran Associates, Inc., Barcelona, 2016, pp. 1019–1027. ISSN 0302-9743. ISBN 9789537619084. doi: 10.1201/9781420049176. http://papers.nips.cc/paper/6241-a-theoretically-grounded-application-of-dropout-in-recurrent-neural-networks.pdf.

37.

Srivastava

Hinton

Krizhevsky

Sutskever

and Salakhutdinov

, Dropout: A Simple Way to Prevent Neural Networks from Overfitting, Journal of Machine Learning Research 15(1) (2014), 1929–1958. ISBN 2013023995. doi: 10.1016/j.micromeso.2003.09.025.

38.

Kingma

D.P.

and Ba

J.L.

, Adam: A Method for Stochastic Optimization, in: 3th International Conference on Learning Representations (ICLR 2015), San Diego, CA, 2015. http://arxiv.org/abs/1412.6980.

39.

Dua

and Graff

, UCI Machine Learning Repository, 2017. http://archive.ics.uci.edu/ml.

40.

Hong

Pinson

Fan

Zareipour

Troccoli

and Hyndman

R.J.

, Probabilistic energy forecasting: Global Energy Forecasting Competition 2014 and beyond, International Journal of Forecasting 32(3) (2016), 896–913. ISBN 0169-2070. doi: 10.1016/j.ijforecast.2016.02.001.

41.

Bergstra

and Bengio

, Random Search for Hyper-Parameter Optimization, Journal of Machine Learning Research 13 (2012), 281–305.

Addressing model uncertainty in probabilistic forecasting using Monte Carlo dropout

Abstract

Keywords

1. Introduction

3. Prediction interval estimation

5.1 Data

Table 5 Average metrics for each dataset and dropout mode (N: No dropout, Y: Regularization dropout, MC: Monte Carlo dropout), for 24 step ahead forecasting

Footnotes

Acknowledgments

References

Table 5
Average metrics for each dataset and dropout mode (N: No dropout, Y: Regularization dropout, MC: Monte Carlo dropout), for 24 step ahead forecasting