While neural networks excel at predicting a single-point estimate of a given target for complex machine learning problems, an open research question is the design of neural generative models able to output a predictive distribution, that can capture the inherent variability of the observations or the model's level of confidence in its predictions. The main motivation behind uncertainty quantification is the design of AI systems for critical applications that are safe, and are mitigating risks while automatizing decision-making. On one hand, Bayesian statistics offer a mathematically grounded framework to reason about uncertainty; however, such models generally require prohibitive computational costs, which make them not widely used in practice. On the other hand, frequentist methods and metrics have been developed for confidence estimation in neural networks, in particular in the classification setting [Brosse et al., 2020] , [Corbière et al., 2019] . Such works address the issue of neural networks calibration [Guo et al., 2017] and detection of out-of-distribution samples [Lee et al., 2018] . But few works focus on generative models based on recurrent neural networks.
On another note, the Transformer model introduced in [Vaswani et al., 2017] has achieved impressive results on sequential data problems. In the field of Natural Language Processing (NLP), Transformers have indeed repeatedly outperformed recurrent neural networks (RNNs), and are now the go-to network architectures to solve complex tasks such as Machine Translation or Language Modeling. The Transformer model is based on a self-attention mechanism, that computes dot-product attention for every element of a sequence with respect to all others to model their dependency. By computing in parallel multiple heads of self-attention and by stacking layers of such multi-head attention, the model learns long-range dependencies better than previous state-ofthe-art models for sequential data.
Since the release of "Attention is all you need" [Vaswani et al., 2017] , the Deep Learning community have been eager to seize the nuts and bolts of the model to improve its training algorithm or to adapt its architecture to new use-cases and types of sequential data, see for instance BERT and its variants ( [Devlin et al., 2019] , [Radford et al., 2018] , [Radford et al., 2019] , [Raffel et al., 2019] ) and [Chen et al., 2018] , [Hao et al., 2019] . Yet, few publications actually focus on the dependency structure at the heart of the Transformer architecture which provides a promising approach to model sequential data. We think that designing a recurrent network inspired from the original Transformer with statistical priors in its architecture could provide a powerful statistical model for capturing the distribution of the observations for sequence prediction problems.
To that end, we introduce the Sequential Monte Carlo (SMC) recurrent Transformer which presupposes that the keys, queries, values and self-attention parameters are unobserved latent states evolving randomly through time. The model relies on a dynamical system inspired from the Transformer, capturing the uncertainty by replacing deterministic self-attention sequences by latent trajectories. The combination of self-attention through time with unobserved model noise allows to generate observations with a complex statistical structure. Self-attention vectors being unobserved stochastic variables, the log-likelihood of the observations is intractable and needs to be estimated. In this paper, we propose to use particle filtering and smoothing methods to draw samples from the distribution of hidden states given observations. The proposed algorithm is based on the auxiliary particle filter of [Liu and Chen, 1998, Pitt and Shephard, 1999] , the most widely used particle filtering method and a generalization of previous approaches such as the algorithms proposed in [Gordon et al., 1993] and [Kitagawa, 1996] . The theoretical properties of such methods to estimate the unknown distributions of the internal states and observations have been widely studied, usually in the context of hidden Markov models, see for instance [Del Moral, 2004] , [Cappé et al., 2005] , [Del Moral et al., 2010] , [Del Moral et al., 2015] , [Douc et al., 2011] , [Dubarry and Le Corff, 2013] , [Olsson et al., 2017] , [Nguyen et al., 2017] .
Fitting the Transformer approach to general state space modeling provides a new promising and interpretable statistical framework for sequential data and recurrent neural networks. From a statistical point of view, the SMC Transformer provides an efficient way of writing each observation as a sophisticated mixture of previous data, while the approximated posterior distribution of the unobserved states captures the states dynamics. We evaluate our model on a series of experiments that first assess the model's ability to predict accurately a known distribution of observations on synthetic datasets, and then apply it to uncertainty quantification in a time-series forecasting problem. The results show that the SMC Transformer manages to capture efficiently the known observation models in the synthetic setting. When performing the task of time-series forecasting on a critical dataset (daily deaths from the Covid-19), the SMC Transformer allows to maintain accurate mean predictions with satisfactory confidence intervals, that take in account the variability of the observations found in each test example.
The Transformer model has been originally developed to propose a new sequence transduction model without recurrence or convolution, and as an alternative to recurrent neural networks for sequence modeling, [Vaswani et al., 2017] . This approach relies entirely on the self-attention mechanism ([Lin et al., 2017] ) to model global dependencies regardless of their distance in input or output sequences. Let (X s ) s 1 be a sequence of observations indexed by N. Transformer models are designed to predict an output X s , for a given index s, from input data X −s . Each input data X s is associated with a query q s and a set of key-value (k s , v s ) computed from linear transformations of the input. From the set of keys and queries, a softmax score function is first computed, which determines how much focus to place on each input in X −s as X s is processed. Transformer models use a scaled dot product attention to compute this score π: π s = softmax(Q s K T s / √ r), where each line of Q s (resp. K s ) contains an input query (resp. keys) and r denotes the dimensionality of keys and queries. Then, this softmax score is used to compute a weighted sum of the values vectors. The final attention is written as:
where V s is the matrix whose rows are the values associated with each input data. The Transformer uses multihead self-attention: the data hidden representations are linearly projected in h subspaces where the attention is computed in parallel, leading to h outputs then concatenated back together.
In real-world machine learning applications, the auxiliary states (queries, keys, values) as well as the observations used to train the neural network are prone to be very noisy. The use of a generative model replacing these deterministic states by random variables evolving according to a dynamical model allows to take into account the uncertainty in the estimation procedure instead of choosing fixed deterministic states. However, considering queries, keys and values as unobserved random variables leads to an intractable likelihood function as the log-likelihood of the observed data X 1:T , where X u1:u2 stands for (X u1 , . . . , X u2 ) for u 1 u 2 , is obtained by integrating out all latent variables which cannot be done analytically. The exact computation of the likelihood function is therefore not possible in general state spaces. Maximum likelihood estimation cannot be performed directly but a gradient descent algorithm may still be defined using Fisher's identity, see for instance [Cappé et al., 2005] :
where θ denotes the unknown parameters of the model, ζ 1:T denotes all the unobserved states, p θ the joint probability distribution of the observations X 1:T and the latent states and E θ the expectation under p θ . In this paper, we propose to estimate the gradient of the log-likelihood of such general state space model using Sequential Monte Carlo methods, i.e. by a set of random samples associated with non negative importance weights. These particle filters and smoothers approximations combine sequential importance sampling steps to recursively update conditional expectations of the form (1) and importance resampling steps to duplicate or discard particles according to their importance weights. By (1), ∇ θ log p θ (X 1:T ) is then approximated by a weighted sample mean of the form
where (ω m n ) 1 m M are nonnegative importance weights such that M m=1 ω m T = 1 and where ξ m 1:T are trajectories approximately sampled from the posterior distribution of ζ 1:T given X 1:T when the parameter is θ. This approximation of the score function can be plugged into any stochastic gradient algorithm to find local minima of θ → − log p θ (X 1:T ).
In the specific case of sequential data, Transformers-based approaches may be introduced to provide a model of the conditional distribution of X t given ∆ past observations X t−∆:t−1 where 1 ∆ t.
In a SMC Transformer with a unique layer, for all 1 s t and all 1 h n heads , define,
Finally, self-attention of the input data is computed, for all 1 s ∆, as,
where π h s (t) denotes the s-th component of π h (t) (i.e. the self attention weight of the observation t−s), (Σ h,z ) 1 h n heads are unknown semi definite-positive matrices and (ε h z ) 1 h n heads are independent standard Gaussian random vectors in R r . Therefore, z h (t) is a Gaussian random variable with mean
The reparametrization trick is used to write z h (t) as a deterministic function of µ h (t) and Σ h with ε h a random variable that does not depend on the parameters. Such trick provides a differentiable transition for the optimization process, see [Kingma and Welling, 2014] . The output r t is then computed with layer normalization and residual connection steps depending on a parameter η state .
In a classification setting, the observation model provides a probability vector G η obs (r t ) on the finite observation space based on the self-attention vectors. In a regression framework, the observation model is given by
where G η obs is a FFNN with linear ouput layer and ε t is a centered noise, for instance a centered Gaussian random vector with unknown variance Σ obs . Let θ be the vector that contains all the unknown parameters of the model:
By Section 3.1, the unobserved state at time t is ζ t = {z(t), q(t), κ(t), v(t)} and the complete-data likelihood may be written
where by convention if t − ∆ 1 then u t−∆:s = u 1:s . The associated probability density functions are
in a classification setting, and
in a regression setting, where ϕ µ,Σ is the Gaussian probability density function with mean µ and covariance matrix Σ. By (1) and (2), the sequential Monte Carlo algorithm approximates ∇ θ log p θ (X 1:T ) by a weighted sample mean:
where the importance weights (ω m T ) 1 m M and the trajectories ξ m 1:T are sampled according to the particle filter described below. In this paper, we propose to estimate all the parameters of the recurrent architecture based on a gradient descent using S M θ,T . All parameters related to the noise (the covariance matrices) are estimated using an explicit Expectation Maximization (EM) update [Dempster et al., 1977] each time a batch of observations is processed, see the appendix materials for all details.
Particle filtering/smoothing algorithm. For all t 1, once the observation X t is available, the weighted particle sample {(ω m t , ξ m 1:t )} N m=1 is transformed into a new weighted particle sample. This update step is carried through in two steps, selection and mutation, using the auxiliary sampler introduced in [Pitt and Shephard, 1999] . New indices and particles {(I m t+1 , ξ m t+1 )} N m=1 are simulated independently as follows:
2. Sample ξ m t+1 using the model introduced in Section 3.1 with the resampled trajectories.
For any m ∈ {1, . . . , N }, the ancestral line ξ 1:t+1 is updated as follows ξ m 1:t+1 = (ξ I m t+1 1:t , ξ m t+1 ) and is associated with the importance weight defined by
Therefore, in a classification setting, ω m t+1 ∝ [G η obs (r m t+1 )] Xt+1 , and in a regression setting, ω m t+1 ∝ exp{− X t+1 − G η obs (r m t+1 ) 2 Σ obs /2}, where for any vector u, u 2 Σ obs = u T Σ −1 obs u. This smoother introduced in [Kitagawa, 1996] (see also [Del Moral, 2004 ] for a discussion) approximates the joint smoothing distributions of the latent states given the observations using the genealogy of the particles produced by the auxiliary particle filter. The genealogical trajectories are defined recursively and updated at each time step with the particles and indices (ξ m k+1 , I m k+1 ). As a result, at each time step, the algorithm selects an ancestral trajectory by choosing its last state at time k which is extended using the newly sampled particle ξ m k+1 . As explained for instance in [Kitagawa, 1996, Kitagawa and Sato, 2001] , [Fearnhead et al., 2010] and [Poyiadjis et al., 2011] , this algorithm suffers from the path degeneracy issue. At each time t 1, the first step to build a new trajectory is to select an ancestral trajectory chosen among M existing trajectories, as the number of resampling steps increases, the number of ancestral trajectories which are likely to be discarded increases. There are many solutions to improve the approximation S M θ,T ; in this paper, we propose to use the fixed-lag smoother of [Olsson et al., 2008] , which means that for each 1 t n, the trajectories ξ m t−∆:t−1 involved in S M θ,T are only resampled up to a few time steps after t.
Inference and predictive distribution. Based on the generative model proposed in this paper, usual objectives are the state estimation problem, which aims at recovering the latent attention parameter z t at time t given the observations X 1:t and the inference problem which aims at approximating the distribution of X t given X 1:t−1 . After a training phase which produces an estimate θ of θ, the state estimation problem is usually solved by approximating the posterior mean of z t given the observations X 1:t when the model is driven by the parameter θ:
which may be approximated using the weighted samples at time t − 1 by
This distribution may be computed using the states dynamics which implies that p θ (z t |ξ m 1:t−1 , X 1:t−1 ) is a Gaussian probability density function. A Monte Carlo approximation of the predictive probability p M θ (X t |X 1:t−1 ) may be obtained straightforwardly by sampling from p θ (z t |ξ m 1:t−1 , X 1:t−1 ). This Monte Carlo estimate can be extended straightforwardly to predictions at future time steps.
Consider first an experimental setting where the observations model is known to assess the ability of our model to capture the distributions of the observations. To that end, we designed two synthetic auto-regressive time-series with a sequence length of 24 observations. For model I, one data sample X = (X 0 , X 1 , ..., X 24 ) is drawn as follows:
where (ε t ) 1 t 24 are i.i.d standard Gaussian variables independent of X 0 . For model II, the law of a new observation given the past is multimodal and drawn as follows: where (ε t ) 1 t 24 are i.i.d standard Gaussian variables independent of X 0 and (U t ) 1 t 24 are i.i.d Bernoulli random variables variables with parameter p independent of X 0 and of (ε t ) 1 t 24 . In the first experiment, the dataset is sampled with α = 0.8 and σ 2 = 0.5. In the second experiment, the dataset is sampled with α = 0.9, β = 0.6α, p = 0.7 and σ 2 = 0.3. More details on the hyperparameters used for such experiments are provided in the appendix. Here, all covariance matrices of the SMC Transformer are assumed to be scalar. Table 1 displays the loss values and their associated standard deviations at the end of training over a 5-fold cross-validation. The displayed losses are the mean squared errors between the predictions (resp. weighted mean of the output particles) and the true observations for the LSTM (resp. the SMC Transformer).
The main interest of our generative model and our estimation procedure is displayed in Figure 1, Table 2 and Figure 2 . For each test example, at each time step t, 1000 samples are drawn from the SMC estimate of the law of X t+1 given X 0:t . A MC-Dropout approach as described in [Gal and Ghahramani, 2015, Section 4] is also used to estimate the uncertainty of X t+1 from 1000 stochastic forward passes through the LSTM network with dropout. In this table, dropout is added after the output of the LSTM layer. Figure 1 illustrates that the samples from the SMC estimate match the true distribution at each time step while the dropout samples highly underestimate its variability. These samples are compared in Figure 1 to the true 95% confidence interval at each time step which is available in this synthetic setting. Our generative model captures the probability distribution of the observations while the LSTM models fail to estimate the known variance of the observations. Based on these 1000 samples, Table 2 provides the empirical estimate of the mean squared error of the predictive distribution of X t+1 given the past for all time steps t. For Model I, it is given by E[(X t+1 − αX t ) 2 |X t ] (and the true value is σ 2 = 0.5) and by
for Model II. The means and standard deviations over the all time steps and all test examples illustrate how the SMC samples capture the true variability of the law of a new observation given the past. In this table, LSTM with dropout also in the LSTM layer are displayed. Dropout networks highly underestimate the variability which yields a predictive MSE too small for both the monomodal distribution (Model I) and the multimodal distribution (Model II). This is confirmed by Figure 2 which shows histograms over the 100 test examples and over all time steps of the frequency of the 1000 samples which fall in the true confidence interval. Over all test samples and time steps, 78.4% (resp. 93.6%) of the SMC with M = 30 (resp. LSTM with d = 64 and dropout = 0.5) samples fall into the 80% confidence interval. Additional results are available in the appendix. 
The performance of the stochastic Transformer is analyzed using Covid-19 data 1 which gathers daily deaths from the Covid-19 disease in 3261 US cities. Cities with less than 100 deaths over the time period considered were discarded from the dataset: all algorithms are thus trained on 886 cities decomposed into 80% for training and 20% for test and validation. Figure 3 displays the 20 days ahead predictions, where the true observation is not available after day 40 and replaced by a sample from the estimated model, for three cities in the test dataset. Here, the variability of the observations is different for every sample: for our model, we first estimated a global variance of the observations at training time, and then fine-tuned the estimated noise per test sample at inference, using 30 iterations of an EM algorithm on the first 40 days. As the time horizon increases, the variability of our predictions increases, thus maintaining a satisfactory combination of mean prediction and error bars: the latter are crucial to model both the uncertainty in the observations measurements (e.g potential errors in data collection), and the natural randomness found in such problem, where daily death rates depend on complex dynamics and a multitude of external factors. On the other hand, the underestimation of the output variability given by the MC-Dropout samples for a rate of 0.2 is highlighted as it leads in some test examples to drifted predictions (see Maryland example on the graph). When increasing the dropout rate up to 0.5, the increase in variability comes at the cost of performance's degradation (see Fairfax example and the training table): this suggests that the confidence interval given by MC-Dropout is ill-calibrated, and that there is no easy tuning of the dropout rate that could give satisfactory error bars. Additional graph, results and full details about the hyper-parameters used for training are available in the appendix. 
The SMC Transformer is part of an emerging line of research which focuses on reconciling traditional model-based approaches and model-free ones that arose from the development of Deep Learning. One such approach similar to the SMC Transformer is the Particle Filter Recurrent Neural Network (PF-RNN) from [Ma et al., 2019] . Our work differs in several ways. First, the importance weights of the PF-RNN do not depend directly on the output of the RNN model, i.e. on an observation model estimated thanks to the SMC approach, but on a exterior learned function. The SMC Transformer uses Fisher's Identity to estimate the gradient of the likelihood of the observations, while the PF-RNN directly optimizes the classic cross-entropy loss and adds an evidence lower bound term in the objective function. The prediction algorithm of the PF-RNN only leverages the SMC algorithm to improve performance at training time and outputs a single-point estimate, thus failing to capture the observations distribution, which is the focus of our paper. On this topic, our work is related to Bayesian Deep Learning: two popular Bayesian methods to quantify uncertainty are MC-Dropout from [Gal and Ghahramani, 2015] and Bayes by Backprop from [Blundell et al., 2015] and its extension to RNNs [Fortunato et al., 2017] . By introducing randomness in the network's parameters or in the training algorithm, they both train an ensemble of neural networks giving a predictive distribution which tends to be overconfident (as illustrated in Section 4); while our method, based on stochastic states, learns directly an observation model almost able to capture perfectly the true variability of the observations.
In this paper, we proposed the SMC Transformer, a novel recurrent network naturally capturing the observations distribution. The model maintains a distribution of self-attention parameters as latent states, estimated by a set of particles. It then outputs a distribution of predictions instead of a single-point estimate, and our inference method gives a flexible framework to quantify the observations variability.
The number of particles of the embedded algorithm is a hyper-parameter that can be tuned depending on the complexity of the task and dataset considered, and by taking in account the trade-off between training time and precision. To our knowledge, this is the first method dedicated to estimate uncertainty in the newly-developed Transformer model, and one of the few focusing on uncertainty quantification in the context of sequence prediction. This SMC Transformer layer could be used as a "plug-and-play" layer for uncertainty quantification in a deeper neural network representing sequential data: the data could be first encoded in a multi-layer neural network before being processed sequentially by the SMC Transformer layer. We specifically chose to focus on the Transformer model by having in mind future applications in the NLP field. We are indeed particularly interested in the diversity in language generation that could arise from the particle filter algorithm at inference. For such task, the distributions of words created by the SMC Transformer could naturally mimic the diversity and richness found in natural language.
The unobserved state at time t is ζ t = {z(t), q(t), κ(t), v(t)} and the complete-data likelihood may be written p θ (X 1:T , ζ 1:T ) = T t=1 p θ (ζ t |ζ t−∆:t−1 , X t−∆:t−1 )p θ (X t |ζ t−∆:t , X t−∆:t−1 ) .
In the regression framework of this paper, the associated probability density functions is
where ϕ µ,Σ is the Gaussian probability density function with mean µ and covariance matrix Σ. A graphical representation of our model which describes the dependency between states and observations is proposed in Figure A .
Hidden layer
Cell input
Output layer Figure 4 : Graphical representation of the SMC transformer for sequential data.
Particle filtering/smoothing algorithm. For all t 1, once the observation X t is available, the weighted particle sample {(ω m t , ξ m 1:t )} N m=1 is transformed into a new weighted particle sample. This update step is carried through in two steps, selection and mutation, using the auxiliary sampler introduced in [Pitt and Shephard, 1999] . New indices and particles {(I m t+1 , ξ m t+1 )} N m=1 are simulated independently as follows:
2. Sample ξ m t+1 using the model with the resampled trajectories.
In the regression framework of the paper, for any m ∈ {1, . . . , N }, the ancestral line ξ 1:t+1 is updated as follows ξ m 1:t+1 = (ξ I m t+1 1:t , ξ m t+1 ) and is associated with the importance weight defined by The algorithm is illustrated in Figure 5 : particles at the last time step are in blue and pink particles are the ones which appear in the genealogy of at least one blue particle. White particle have not been selected to give birth to a path up to the last time. In Figure 5 , the N = 3 genealogical trajectories are ξ 1 0:4 = (ξ 3 0 , ξ 2 1 , ξ 2 2 , ξ 3 3 , ξ 1 4 ), ξ 2 0:4 = (ξ 3 0 , ξ 1 1 , ξ 3 2 , ξ 2 3 , ξ 2 4 ), ξ 3 0:4 = (ξ 3 0 , ξ 2 1 , ξ 2 2 , ξ 3 3 , ξ 3 4 ).
The training algorithm The sequential Monte Carlo algorithm approximates ∇ θ log p θ (X 1:T ) by a weighted sample mean:
where the importance weights (ω m T ) 1 m M and the trajectories ξ m 1:T are sampled according to the particle filter described below. Thanks to Fisher's identity, this approximation only requires to compute the gradient of the state model loglikelihood θ → log p θ (ξ m t |ξ m t−∆:t−1 , X t−∆:t−1 ) and the gradient of the observation model loglikelihood θ → log p θ (X t |ξ m t−∆:t , X t−∆:t−1 ). There is no need to compute the gradient of the weights ω m T which depend on the parameter θ. Using TensorFlow function tf.stop gradient on these weights, this allows to train the model with the following loss function:
log p θ (ξ m t |ξ m t−∆:t−1 , X t−∆:t−1 ) + log p θ (X t |ξ m t−∆:t , X t−∆:t−1 ) .
In this paper, we propose to estimate all the parameters of the recurrent architecture based on a gradient descent using S M θ,T . All parameters related to the noise (the covariance matrices) are estimated using an explicit Expectation Maximization (EM) update [Dempster et al., 1977] each time a batch of observations is processed. For each sequence of observations, the EM update relies on the approximation of the intermediate quantity
log p θ (ζ t |ζ t−∆:t−1 , X t−∆:t−1 ) + log p θ (X t |ζ t−∆:t , X t−∆:t−1 )|X 1:T ] by the following particle-based estimator:
log p θ (ξ m t |ξ m t−∆:t−1 , X t−∆:t−1 ) + log p θ (X t |ξ m t−∆:t , X t−∆:t−1 ) .
Then, Q M θ,T may be maximized with respect to all covariances to obtain the new estimates. This is a straightforward update which yields for instance for Σ obs for the p-th update:
where r m t are the resampled particles at time t. The estimate Σ obs of Σ obs is then updated as follows:
where η p is a learning rate chosen by the user (η p = p −0.6 in the experiments).
Inference and predictive distribution used in the experiments. To solve the inference problem, note that p θ (X t |X 1:t−1 ) = p θ (X t , z 1:t |X 1:t−1 )dz 1:t ,
which may be approximated using the weighted samples at time t − 1 by
A Monte Carlo approximation of the integral involved in the predictive probability p M θ (X t |X 1:t−1 ) may be obtained straightforwardly. In the experiments, p M θ was approximated as follows. For all 1 k N (with N = 1000 in our experiments), sample independently I k in {1, . . . , M } with weights proportional to (ω m t−1 ) 1 m M and then z k t from p θ (z t |ξ I k 1:t−1 , X 1:t−1 ), i.e. from our transformer model with the selected past trajectory of hidden states ξ I k 1:t−1 . Finally r k t can be computed from z k t and p M θ is approximated by N −1 N k=1 ϕ Gη obs (r k t ), Σ obs (X t ) i.e. by a mixture of Gaussian distributions with means given by (r k t ) 1 k N and variance Σ obs .
The synthetic datasets were generated with 1000 samples: we used 800 of them for training and 100 of them for test and validation. For training the models, we used a batch size of 32, a number of epochs equal to 50, and the ADAM algorithm with a learning rate of 0.001 and the original custom schedule found in [Vaswani et al., 2017] . Table 4 displays the loss values and their associated standard deviations at the end of training over a 5-fold cross-validation. The displayed losses are the mean squared errors between the predictions (resp. weighted mean of the output particles) and the true observations for the LSTM (resp. the SMC Transformer). Additional dropout rates and a Transformer with d = 8 were considered to complete the experiments given in the main part of the paper. These neural networks were also considered in Table 5 .
The hyper-parameters used for training the SMC Transformer and the LSTM with dropout were the following: the models were trained during 50 epochs with a batch size of 32, with a learning rate of 0.001 for the LSTM and the original learning rate with custom schedule from [Vaswani et al., 2017] . The LSTM with dropout include a dropout layer inside the LSTM layer, and before the output layer. Figure 6 displays another example of 20 days ahead predictions, where the true observation is not available after day 40 and replaced by a sample from the estimated model, for two cities in the test dataset. The figure shows that the confidence intervals given by the SMC Transformer capture almost entirely the true sequence of predictions for these two cities with high variability in daily death rates, which is not the case for the two LSTM with dropout. Figure 6 : 20 days ahead predictions of covid-19 daily deaths for 2 cities. The black lines represent the true predictions, and the colored areas represent the confidence interval given by a SMC Transformer with 10 particles, while the grey stripped (resp. dotted) areas represent the confidence intervals given by a LSTM with a dropout rate equal to 0.2 (resp. to 0.5). 