arXiv CSopen access

Unveiling Stochasticity: Universal Multi-modal Probabilistic Modeling for Traffic Forecasting

Unveiling Stochasticity: Universal Multi-modal Probabilistic Modeling for Traffic Forecasting Weijiang Xionga , Robert Fonoda , Nikolas Geroliminisa,∗

arXiv:2604.16084v1 [cs.LG] 17 Apr 2026

a

Urban Transport Systems Laboratory (LUTS), EPFL, Lausanne, 1015, Switzerland

Abstract Traffic forecasting is a challenging spatio-temporal modeling task and a critical component of urban transportation management. Current studies mainly focus on deterministic predictions, with limited considerations on the uncertainty and stochasticity in traffic dynamics. Therefore, this paper proposes an elegant yet universal approach that transforms existing models into probabilistic predictors by replacing only the final output layer with a novel Gaussian Mixture Model (GMM) layer. The modified model requires no changes to the training pipeline and can be trained using only the Negative LogLikelihood (NLL) loss, without any auxiliary or regularization terms. Experiments on multiple traffic datasets show that our approach generalizes from classic to modern model architectures while preserving deterministic performance. Furthermore, we propose a systematic evaluation procedure based on cumulative distributions and confidence intervals, and demonstrate that our approach is considerably more accurate and informative than unimodal or deterministic baselines. Finally, a more detailed study on a real-world dense urban traffic network is presented to examine the impact of data quality on uncertainty quantification and to show the robustness of our approach under imperfect data conditions. Code available at [link provided upon acceptance]. Keywords: Traffic forecasting, Uncertainty quantification, Deep learning, Gaussian mixture models ∗

Corresponding author. Email addresses: [email protected] (Weijiang Xiong), [email protected] (Robert Fonod), [email protected] (Nikolas Geroliminis)

1. Introduction Spatio-temporal modeling is a fundamental task in machine learning with applications in various domains, such as sustainable power systems (Bessa et al., 2014), climate change (Maclean, 2020), transportation systems (Yuan and Li, 2021) and video data analysis (Zhou et al., 2025). Specifically, traffic forecasting has been an active research topic in this domain, owing to its vital role in smart cities (Tedjopurnomo et al., 2022; Chen and Chen, 2019). With a traffic forecasting model, future traffic states of the city can be predicted using historical data, providing essential grounds for adaptive traffic management (Ravish and Swamy, 2021; Wang et al., 2018). In earlier research, statistical methods, such as Historical Average (HA) and Autoregressive Integrated Moving Average (ARIMA), were developed as baseline predictors that use the periodic patterns and temporal dynamics of traffic (Vlahogianni et al., 2014). However, these approaches have limited performance because traffic data contains complex dependencies in both spatial and temporal dimensions: the future state of one road depends on its own dynamics and on nearby roads (Ermagun and Levinson, 2018). As a result, the modern research trend focuses on learning these correlations to improve prediction accuracy (Jiang and Luo, 2022). Since the collected data is deterministic, i.e., only one value is measured for a traffic state variable at a specific location and time, the most natural evaluation is to compare error metrics, such as Mean Absolute Error (MAE) (Ferreira et al., 2023). These facts and practices have shaped the paradigm for developing and evaluating traffic forecasting methods, and traffic forecasting has long been framed in a deterministic setting, i.e., one future value is predicted and evaluated. However, urban traffic is highly dynamic and uncertain due to various factors, such as traffic signals, pedestrian activities, weather conditions, traffic incidents, and the randomness in drivers’ behaviors (Mauro et al., 2015). As a result, the process that generates deterministic traffic data is stochastic, and the collected data can be regarded as one realization of that process. From the perspective of modeling, the future traffic variables should be modeled with probabilistic distributions instead of deterministic values (Ghosh et al., 2010). On the other hand, real-world decision-making processes, such as traffic control, should consider all possibilities to provide safe and robust solutions covering all outcomes (Wanke et al., 2004). Moreover, the randomness in traffic can have a big impact and result in distinct outcomes. For example, travel time on an arterial road can have multiple probable values 2



'DWD *DXVVLDQ 0XOWL0RGDO 'HWHUPLQLVWLF ,QWHUPHGLDWH6SHHGV

,QWHUPHGLDWHVSHHGV KDYHORZHUGHQVLW\

3UREDELOLW\'HQVLW\

    









7UDIILF6SHHGNPK





Figure 1: Multi-modal data distribution and various predictions

due to uncertain waiting times at traffic signals (Ramezani and Geroliminis, 2012). Consequently, the data distribution can be multi-modal, i.e., having multiple peaks in the probability density function. Figure 1 illustrates different predictive methods for a multi-modal traffic speed distribution, where a road segment can be either congested or freeflow and intermediate states are less common. This scenario is realistic when the traffic is highly uncertain (especially for urban traffic) or when existing information is insufficient to distinguish the two states (e.g., for long-term predictions). In such a case, predicting a deterministic estimate can conceal the uncertainty and lead to false confidence, even though it can have low deterministic error metrics when averaged over multiple trials. A unimodal distribution (e.g., Gaussian) can indicate the level of uncertainty with a large variance; however, this prediction can be too vague to provide useful information. Building on these insights, this paper proposes to model traffic variables as multi-modal random variables. We create a universal and effective approach to enhance traffic forecasting models by expanding the prediction scope from a most-likely value to multi-modal probabilistic distributions. We show that the adaptation requires only a minimal modification: replacing the last output layer with our novel Gaussian Mixture Model (GMM) layer, yet the performance gain and robustness are consistently demonstrated through experiments on various datasets. With the expressive power of multi-modal 3

predictive distributions, this approach can unveil the underlying stochasticity in deterministic data and provide a stronger basis for decision-making processes that rely on these predictions, e.g., traffic signal control in transportation systems (Tsitsokas et al., 2023). Despite the simplicity of the model modification, improvements in probabilistic evaluation metrics are significant, without sacrificing deterministic prediction performance. In summary, the main contributions of this paper are: • A simple yet effective approach to enhance the descriptive power of traffic forecasting models by expanding the prediction space into multimodal probabilistic distributions. • A general and comprehensive evaluation framework for multi-modal probabilistic predictors, applicable whenever probability density is available. • Extensive experiments on multiple traffic forecasting datasets to demonstrate the advantages of multi-modal probabilistic predictions. • A case study to show the impact of data quality on uncertainty and highlight the robustness of our approach. 2. Related Work This section provides a brief review of the two most related research topics, i.e., traffic forecasting and uncertainty quantification, since our approach lies at their intersection. 2.1. Traffic Forecasting Traffic forecasting has been a widely studied research topic due to its significance both as a complex spatio-temporal modeling problem and as a practical application in transportation systems (Vlahogianni et al., 2014). Most modern approaches are based on deep learning, where the neural network models are designed to cope with the dependencies in both spatial and temporal dimensions (Jiang and Luo, 2022). The Diffusion Convolutional Recurrent Neural Network (DCRNN) is a pioneering work that learns the spatial relations through a diffusion convolution over the graph structure of the transportation network and captures the temporal correlations using a Recurrent Neural Network (RNN) (Li et al., 2018). To improve the efficiency 4

of temporal modeling, the Spatio-Temporal Graph Convolutional Networks (STGCN) applies a convolution over the time dimension, thus avoiding the sequential nature of RNNs (Yu et al., 2018). In fact, traffic conditions at different locations are often intertwined, and data correlations often go beyond geographical distances (Kwak et al., 2021). Therefore, many follow-up works have focused on adaptively learning relations from data, e.g., by incorporating various attributes (Xiao et al., 2026) or by spatio-temporal attention (Liu et al., 2025). The Graph WaveNet (GWNet) (Wu et al., 2019) proposes enhancing the pre-defined graph with an adaptive graph, which is constructed by the pairwise product of learnable node embeddings. The Attention-based Spatial-Temporal Graph Convolutional Networks (ASTGCN) (Guo et al., 2019) further introduces a more dynamic attention mechanism (Vaswani et al., 2017) for spatial and temporal dimensions, which effectively assigns weights to the spatial or temporal neighbors according to the learned correlation patterns. Similarly, the Graph Multi-Attention Network (GMAN) (Zheng et al., 2020) employs an encoderdecoder structure with an attention mechanism on both spatial and temporal dimensions. In addition, the encoder messages are enhanced by a transform attention between past and future time steps, allowing the model to alleviate error propagation. The feature embedding technique also plays an important role in traffic forecasting. The Spatial and Temporal Identity Information (STID) (Shao et al., 2022) identifies that a bottleneck in modern traffic forecasting models is the lack of sufficient information to distinguish different samples, and proposes augmenting feature embeddings with metadata, such as time of day and day of week. The Multivariate Time Series Graph Neural Network (MTGNN) (Wu et al., 2020) designs a graph learning module that combines the node embeddings with static features to formulate more comprehensive node descriptions. Additionally, MTGNN organizes message exchanges for neighbors at different distances, enabling larger receptive fields compared to the traditional GNNs. A similar idea has also been utilized in the Dynamic Multihop Graph Attention Network (DMGAN) (Li et al., 2022). Later, a SpatioTemporal Adaptive Embedding (STAE) technique has been introduced to augment the traffic value embeddings with a set of learnable adaptive embeddings (Liu et al., 2023). As a result, the STAE technique can significantly improve the performance of a vanilla Transformer (Vaswani et al., 2017).

5

2.2. Uncertainty Quantification Research on uncertainty quantification falls within the context of Bayesian Deep Learning, where uncertainty is generally categorized as aleatoric and epistemic uncertainties (Kendall and Gal, 2017). Epistemic uncertainty reflects insufficient knowledge to determine the model parameters, and can be reduced by a larger collection of training data. In contrast, aleatoric uncertainty refers to the inherent stochasticity in data, which originates from the underlying physical data generation processes. Therefore, these two uncertainties are often modeled by probabilistic distributions over the model parameters and output space respectively. Earlier studies investigated uncertainty in arterial road travel time distributions (Ramezani and Geroliminis, 2012), or the multi-regime nature of traffic speed with Mixture-of-Expert formulations (Coric et al., 2011). In addition, the idea of modeling uncertainty has also been studied in other time series forecasting tasks outside transportation domain, such as wind power prediction (Zheng and Zhang, 2023) and medical data analysis (Tan et al., 2020). With the rise of deep learning, Bayesian ideas such as Monte Carlo Dropout (Gal and Ghahramani, 2016) were combined with traffic forecasting models, for example, placing distributions (Wu and James, 2021) over the weights of a STGCN (Yu et al., 2018). Similar practices have also been proposed to develop a Bayesian Graph Convolution Network for spatial uncertainty modeling (Fu et al., 2021). Later, both epistemic (weight) uncertainty and aleatoric uncertainty have been combined in one traffic prediction model (Qian et al., 2024), and then unified together with adaptive weight averaging (Qian et al., 2023). Frequentist approaches have emerged in parallel, for example, the quantile regression formulation has been applied to traffic prediction to adapt existing models into quantile predictors (Wu et al., 2023). To improve the generalization of uncertainty prediction, a spectral-normalization technique has been proposed to stabilize the training of Bayesian RNNs (Sengupta et al., 2024). Traffic forecasting and imputation have also been combined under one joint framework, where a copula-based multivariate distribution is learned to handle the joint distribution of multiple traffic variables (An et al., 2025). The main limitation of most existing methods is that they restrict prediction to a deterministic value or a unimodal (e.g., Gaussian) distribution, which fails to capture the complex stochasticity of traffic. On the contrary, 6

our approach features the capability to explicitly predict multi-modal distributions, which can express multiple possible but distinct future traffic states. Our modular solution can be easily integrated into various model architectures while keeping the training pipelines identical. Additionally, our evaluation framework provides comprehensive metrics to evaluate the whole predictive distribution and confidence intervals at multiple confidence levels. 3. Methodology This section presents our straightforward approach to transform deterministic traffic forecasting models into effective probabilistic predictors using GMMs, and introduces the evaluation metrics to assess the quality of predicted distributions. 3.1. Multi-Modal Probabilistic Traffic Forecasting The goal of traffic forecasting is to predict the future traffic states at multiple future time steps over the transportation network, given the historical traffic states. Since the transportation network naturally introduces spatial dependencies among locations, the traffic forecasting problem has often been modeled as a spatio-temporal prediction task. In general, the structure of the neural networks can be described as follows: Z = ST_Backbone (Xt−Th :t , G) Xt+1:t+Tf = Det_Pred (Z) .

(1)

ST_Backbone(·) is a spatio-temporal modeling neural network (e.g., DCRNN (Li et al., 2018) or GraphWaveNet (Wu et al., 2019)) that extracts hidden representations Z ∈ RN ×D using historical data and graph structure, where D is the feature dimension. In the input to the backbone, Xt−Th :t ∈ RN ×Th represents the traffic states for all N locations across Th historical time steps up to the current time t. Similarly, Xt+1:t+Tf ∈ RN ×Tf represents the future traffic states for all locations up to Tf steps ahead of the current time t. Following Li et al. (2018), Tf and Th are often set to an equal value T in practice. G = (V, E) is the graph structure of the transportation network, whose nodes V = {v1 , v2 , . . . , vN } are the locations with collected traffic data, and edges E indicate the spatial dependencies among the locations. The edges can often be equivalently transformed into an adjacency matrix A ∈ RN ×N , which can be binary for unweighted graphs and real-valued for weighted graphs. 7

Det_Pred(·) is a deterministic layer that predicts the future traffic state values using the hidden features Z from the backbone. Most forecasting methods predict T future time steps in one forward pass (Lim and Zohren, 2021), which is free from error accumulation in autoregressive multistep predictions. As a result, the Det_Pred layer can be described as: Xt+1:t+Tf = Linear(Z) = WZ⊤ + b,

(2)

where W ∈ RT ×D and b ∈ RT are the learnable weights and biases independently applied to each location. As discussed in Section 1, a major limitation of the deterministic paradigm is that the model is always trained to predict a single most-likely value even when there are multiple possibilities. Concretely, a road segment may develop into a congested case with low speed or a free-flow case with high traffic speed, which cannot be determined with current information. Restricting the prediction space to a single value prevents the model from describing such complex uncertainty and hinders its ability to capture the stochastic nature of the data. This limitation in descriptive power also creates challenges for prediction-based decision-making, as deterministic predictions may create an illusion of certainty and undermine system reliability. Therefore, this paper proposes to capture the diverse possibilities by modeling the traffic forecasting problem from a probabilistic perspective: Z = ST_Backbone (Xt−Th :t , G) p(Xt+1:t+Tf ) = GMM_Pred (Z) ,

(3)

where p(Xt+1:t+Tf ) is the predicted probability distribution (instead of a point estimation) of future traffic states. Despite the fundamental change in the prediction space, the adaptation can be achieved in a surprisingly straightforward approach. As illustrated in Figure 2, adapting a deterministic model requires no change in the spatio-temporal backbone. The only modification is to replace the last output layer for deterministic prediction (Det_Pred) with a GMM layer for multi-modal probabilistic predictions (GMM_Pred) 1 . The GMM_Pred layer is a lightweight probabilistic enhancement module designed to be a universal fit for any existing backbone network. Compared to the deterministic method, the GMM layer predicts a continuous 1

A Gaussian Mixture Model is called multi-modal because its probability density function can have multiple peaks, representing multiple modes.

8

𝑝

𝑣

𝑣 GMM Dist.

GMM _Pred

Spatio-Temporal Modeling

Spatio-Temporal Modeling

Var

Det_Pred

Mean

𝑡 Mixing

𝑡

Projection

Spatio-Temporal Input Data

Figure 2: Transforming a deterministic model (left) to a probabilistic predictor (middle) with a GMM prediction layer (right).

multi-modal distribution, endowing the model with the ability to represent not only uncertainty around a mean value but also multiple distinct modes. For simplicity, we assume element-wise independent marginal distributions, which means the prediction at each location and time step is an independent one-dimensional Gaussian mixture distribution: p(x) =

K X

πk N (x; µk , σk2 ),

(4)

k=1

where K is the number of components, N is the Normal density function, and πk , µk , and σk2 are the predicted mixing coefficient, mean, and variance of the k-th component, respectively. Notably, the mixing coefficients πk must sum to 1 to ensure the mixture is a valid probability distribution. For the training loss, we use the GMM Negative Log-Likelihood (NLL)

9

of future values (i.e., labels), whose element-wise form is: − log p(x) = − log = − log

K X k=1 K X k=1 K X

πk N (x; µk , σk2 )

exp log πk N (x; µk , σk2 )

(5)

 (x − µk )2 1 2 − log(2πσk ) , = − log exp log πk − 2 2σ 2 k k=1 

where x is the value of an observed traffic variable at a specific location and time step. The element-wise loss is then averaged over all locations and time steps to obtain the final training loss, without any auxiliary terms. 3.2. Details of the GMM Layer As illustrated in Figure 2, the GMM layer starts with a projection layer to transform location-specific features, implemented as a Linear layer: Z′ = Linear(Z),

(6)

where a linear transform is independently applied to each location, and Z′ ∈ ′ RN ×D is the transformed feature matrix with hidden dimension D′ . Then, Z′ is fed into three separate branches for the mixing coefficients, means, and variances to build the mixture distribution. π = Softmax (Reshape (Linear(Z′ ))) , µ = Reshape (Linear(Z′ )) · s + r,
log σ 2 = Reshape (Linear(Z′ )) .

(7)

The Linear(·) layer in each branch predicts the corresponding parameters for all K components at all T future time steps simultaneously, resulting in output shapes of RN ×(T ×K) for each branch. Then, the Reshape(·) operation reshapes the output into RN ×T ×K . To ensure the mixture is a valid probability distribution, the mixing coefficients π ∈ RN ×T ×K are normalized over the K components using the Softmax function, whose element-wise formula (at a specific location and time step) is: exp(zk ) , (8) π k = PK exp(z ) j j=1 10

Data 3-𝜎 Range 𝜇! range 𝜇! range 𝜇" range -3

-3+s -3+2s

-3+ks

3

Figure 3: Mean value formulation and data range

where zk is the output of the linear layer in the mixing-coefficient branch for the k-th Gaussian component. Following Kendall and Gal (2017), the GMM layer predicts log-variance instead of the variance itself for better numerical stability, i.e., log (σ 2 ) ∈ RN ×T ×K . The mean values µ ∈ RN ×T ×K are reparameterized by the predicted relative offsets together with pre-defined (not learnable) reference values r ∈ RK and scale s ∈ R. The · between scale s and transformed features is a simple multiplication that applies the same scale factor to all locations and time steps. This formulation can inject prior knowledge about the data distribution and stabilize training. In our implementation, both the inputs and labels in Equation 5 are normalized by a Z-score transformation to have zero mean and unit variance, and the training loss is computed in normalized space. The final form of the GMM NLL loss utilizes the log-sum-exp approach to avoid numerical instability. In addition, the GMM layer is initialized to predict a weakly informative prior (Gelman et al., 1995), where all mixing coefficients are equal, the means are the reference values (zero offsets), and the variances are set to one. This initialization involves both neural network parameters and the reference values and scale in the mean-value formulation. Specifically, the linear layer of the mixing branch is initialized with zero weights and bias 1/K, while the mean and variance branches are both initialized with zero weights and bias. Figure 3 illustrates the intuition for the mean value formulation. The reference values r and scale s are chosen so that the predicted mean µ of the farthest component can reach the 3-σ range of the normalized data space more easily. We assume a unit-range output space for the linear layer (e.g., [−1, 1]), and overlap half of the scaled range (e.g., [−s, s]) for two adjacent components (e.g., range of µ1 and µ2 ). Since the data is normalized, the 3-σ 11

range is approximately [−3, 3], and with this formulation we have: −3 + (K + 1)s = 3 → (K + 1)s = 6.

(9)

We use a default of K = 5 mixture components in the GMM layer to allow more flexibility in modeling the data distribution, and therefore s = 1. Then, the corresponding reference values are set to be r = [−2, −1, 0, 1, 2]. 3.3. Derived Predictions From the predicted distribution, we extract point estimates and confidence intervals to demonstrate the versatility of our multi-modal probabilistic predictions across various application scenarios. The Bayesian average of the predicted distribution can serve as a point estimate, which is the weighted sum of the component means: x̂ =

K X

π k µk .

(10)

k=1

We propose a computational method to derive confidence intervals that prioritize high-density regions. The procedure is described in Algorithm 1. Figure 4 illustrates an example confidence interval using a three-component GMM. For intuition, one can imagine a density threshold, i.e., the horizontal line in Figure 4, intersecting the density curve. This threshold line may have multiple intersections with the density curve, and each resulting segment corresponds to a sub-interval at the given confidence level. When the Gaussian components are well-separated and have small variances, the density curve for a GMM with K components can have K peaks. As a result, the threshold line can have up to K line segments inside the density curve, yielding up to K sub-intervals for a given confidence level. 3.4. Evaluating Probabilistic Predictions Since a probabilistic prediction contains much more information than a deterministic counterpart, ideal scoring rules should evaluate the entire distribution rather than derived properties. For this purpose, we follow Tang and Matteson (2021) and use the Continuous Ranked Probability Score (CRPS) (Gneiting and Raftery, 2007) for its joint consideration of calibration, resolution, and uncertainty (Candille and Talagrand, 2005). By definition, the CRPS score evaluates the difference between a predicted Cumulative 12

Probability Density

0.25

Mixture Density Component 1 Component 2 Component 3 Density Threshold 70%_interval

0.20 0.15 0.10 0.05 0.00

4

2

0 Data Value

2

4

Figure 4: Example GMM Confidence interval

Distribution Function (CDF) F and an observed value y: Z ∞ (F (x) − H(x − y))2 dx, CRPS(F, y) =

(11)

−∞

where H is the Heaviside step function, i.e., a step function jumping from 0 to 1 at the label value y. We interpret a deterministic prediction as a Dirac Delta function. Therefore, F (x) reduces to a step CDF at the predicted value, and the CRPS score becomes the absolute error between the prediction and the label. In addition to CRPS, we introduce two metrics to specifically evaluate the predicted confidence intervals, i.e., the mean Average Width (mAW) and the mean Confidence Calibration Error (mCCE). A properly predicted confidence interval should be informative. Therefore, we measure the width of the predicted interval (possibly with multiple sub-intervals indexed by i) at a confidence level, and average it over all locations and time steps to obtain an Average Width (AW). ! 1 X X i i |utv − ltv | , (12) AW = |T ||V | t,v i where |T | and |V | are the number of prediction time steps and locations, and i uitv and ltv are the upper and lower bounds of the i-th predicted sub-interval at time step t and location v. Then, AW is averaged over all confidence levels

13

Algorithm 1 Confidence Interval (CI) Derivation 1: Input: GMM {πk , µk , σk }K k=1 , confidence level c 2: Output: Lower and upper bounds of CI(s) {lk , uk }k 3: Calculate numerical density p(x) at multiple points x.

x has step size ∆x and covers the range of interest

4: Sort p(x) by descending order of probability density 5: Compute a cumulative sum for the sorted density

F (x) = cumsum(p(x) ∗ ∆x) 6: Normalize: F (x) ← F (x)/ max(F (x)) 7: Search the insertion index i for c in F (x) 8: All elements to the left of i are selected (high-density) 9: Backtrack to the original p(x) to find selected x values 10: Find {lk } and {uk } as follows: {lk }: selected x’s where x − ∆x is not selected {uk }: selected x’s where x + ∆x is not selected 11: Sort {lk } and {uk } and pair by index 12: Return {lk , uk }k of interest to obtain mAW 2 . mAW =

1 X AWc , |C| c

(13)

where |C| is the number of confidence levels, and AWc is the average width at confidence level c. While the mAW measures how informative the confidence intervals are, the predicted intervals should also be well-calibrated, e.g., a 90% confidence interval should contain the true label 90% of the time. When the confidence and coverage are misaligned (either over-confident or under-confident), the predicted intervals can mislead decision processes relying on those predictions. In fact, a confidence interval can be regarded as a binary classification, where the model predicts that the label value is within the interval with a certain confidence level. Therefore, we draw inspiration from Guo et al. (2017) and compute the Confidence Calibration Error (CCE) at a confidence level as the absolute difference between the expected confidence level and the 2

We consider 10 confidence levels from 0.5 to 0.95 with a step size of 0.05

14

empirical coverage: ! 1 X X i CCE = 1(ytv ∈ [ltv , uitv ]) − c , |T ||V | t,v i

(14)

where 1(·) is the indicator function that returns 1 if the condition (ytv ∈ i [ltv , uitv ]) is satisfied, and 0 otherwise. In case of multiple sub-intervals, the true label is considered to be within the interval if it falls into any of the sub-intervals. Similar to the metric mAW, we also average the CCE over the same set of confidence levels to obtain the mCCE: 1 X CCEc . (15) mCCE = |C| c 4. Experiments This section first presents the experiment setup, introduces the evaluation results using probabilistic metrics and traditional deterministic metrics, and finally discusses the impact of data quality on uncertainty through an additional study. 4.1. Datasets and Experiment Setup The experiments and evaluations are conducted on two widely recognized datasets, METR-LA and PEMS-Bay (Li et al., 2018), and a novel dataset SimBarcaSpd, derived from the SimBarca dataset with dense urban traffic from 1570 road segments in the 15 km2 area of central Barcelona (Xiong et al., 2025). The dense connections create more complex spatial dependencies, and thus making the urban traffic forecasting much more challenging than traditional datasets with sparser road networks. Contrary to METR-LA and PEMS-Bay, which are collected from highway loop detectors, SimBarcaSpd is derived from a microscopic traffic simulator using the dense road network of central Barcelona. All three datasets provide traffic speed data, whose maximum values differ due to regional speed limits (the minimum is naturally 0). At the same time, error metrics are computed in the same units as those maximum values. These specifications are summarized in Table 1. To show that our approach can be applied to a wide range of traffic forecasting models, we conduct experiments on several representative models, including GWNet (Wu et al., 2019), MTGNN (Wu et al., 2020), STID (Shao 15

Table 1: Dataset specifications

Dataset METR-LA PEMS-Bay SimBarcaSpd

In-Step Out-Step Step Size #Nodes 12 12 10

12 12 10

5 min 5 min 3 min

207 325 1570

Max

Area

70 mph 477 km2 85 mph 325 km2 14 m/s 15 km2

et al., 2022), STAEformer (Liu et al., 2023), and a simple baseline model LSTM-GCN (denoted as LGC), as shown in Figure 5. The LGC model first encodes the temporal dimension using a three-layer LSTM (Hochreiter and Schmidhuber, 1997) and then processes spatial dependencies using a threelayer Graph Neural Network implemented with GCNConv (Kipf and Welling, 2016). Then, the spatial and temporal features are concatenated and passed to a linear layer. GNN

Data

LSTM

Linear

Pred.

Figure 5: Structure of the LGC model

For each of these models, three variants are considered: a deterministic variant (Det), a Normal variant (Norm) predicting a Gaussian distribution, and a GMM variant with the adaptation method in Figure 2. In fact, the Norm variant is a special case of the GMM variant with a single component. For a fair comparison, we utilize a consistent training and evaluation pipeline across all models and variants. The models are trained for 50 epochs with batch size 32, except STAEformer, which is trained for 20 epochs with batch size 8 to fit memory constraints and prevent heavy overfitting. The AdamW optimizer (Loshchilov and Hutter, 2019) is used with a learning rate of 0.0005, weight decay of 0.0001 and beta parameters of (0.9, 0.999). The learning rate gradually increases from 0 to 0.0005 in the first 2 epochs, and decays to 10% and 1% of the chosen learning rate at 75% and 85% progress of total training epochs. In evaluation, the numerical probability density in Algorithm 1 is computed at 500 evenly spaced points from 0 to the maximum 16

speed of each dataset in Table 1. All experiments are conducted on a single Nvidia L40S GPU with 48GB memory, and the implementation is based on PyTorch. 4.2. Overall Evaluation Results Table 2 presents the probabilistic performance metrics, i.e., CRPS, mAW and mCCE, for all three datasets. The metrics are averaged over all prediction horizons to provide a joint evaluation of model performance at both shorter and longer horizons. The rows are grouped by the architecture of base models, with individual variants (Det, Norm and GMM) listed separately. As discussed in the methodology, the CRPS for the Det variant is the same as the Mean Absolute Error (MAE) in the deterministic prediction task. Moreover, mAW and mCCE for the Det variants are not applicable, as these variants do not provide confidence intervals. On the primary metric, CRPS, a consistent hierarchy is observed across all datasets and all base model architectures, where the GMM variant achieves the best performance, followed by the Norm variant, and the Det variant has the worst performance. The difference between the tiers is significant, which shows the advantage of probabilistic predictions over deterministic ones, and the benefit of multi-modal predictions. STAEformer (briefly noted as STAE) achieves the best CRPS scores on all three datasets, due to its advanced spatio-temporal modeling architecture. Another notable observation is that the performance gap between different variants is more pronounced than the gap induced by base model architectures. For example, the variant STAEformer-GMM achieves a CRPS of 2.23 on METR-LA, which is much better than its Norm variant (2.37) and Det variant (3.04). Still, the GMM variant of simpler models, such as LGC and GWNet, can also achieve similar CRPS scores at 2.25 and 2.26, respectively. This observation shows that the probabilistic modeling approach has a stronger impact on prediction performance. The experiments clearly demonstrate the effectiveness and generality of our proposed GMM adaptation method across different backbone architectures. Using the deterministic variant as a baseline (100%), we observe in Table 3 a consistent and significant average relative CRPS improvement of 27.6% by introducing the GMM layer. The confidence intervals are evaluated using both mAW and mCCE, and ideal predictions should have low values for both metrics simultaneously. On METR-LA and SimBarcaSpd, the GMM variants of all models achieve 17

Table 2: Probabilistic performance metrics. Lower values are preferred for all, best result for a backbone in bold.

mCCE

mAW

CRPS

mCCE

mAW

CRPS

SimBarcaSpd

Det 3.11 - 1.65 - 1.05 Norm 2.49 9.53 0.08 1.29 5.06 0.01 0.81 3.07 0.04 GMM 2.25 7.04 0.01 1.22 3.73 0.09 0.75 2.49 0.01 Det 3.31 - 1.67 - 1.06 Norm 2.68 10.58 0.10 1.31 5.74 0.05 0.82 3.29 0.06 GMM 2.38 7.53 0.02 1.22 4.14 0.05 0.76 2.59 0.01

STAE MTGNN GWNet

LGC

mCCE

PEMS-Bay

STID

mAW

METR-LA CRPS

Method

Det 3.13 - 1.58 - 1.04 Norm 2.52 9.67 0.08 1.23 5.09 0.02 0.81 3.18 0.05 GMM 2.26 6.99 0.01 1.14 3.75 0.07 0.76 2.63 0.00 Det 3.12 - 1.57 - 1.00 Norm 2.54 9.92 0.09 1.22 5.17 0.02 0.78 3.12 0.06 GMM 2.27 7.11 0.01 1.15 3.76 0.07 0.72 2.47 0.01 Det 3.04 - 1.56 - 0.98 Norm 2.37 8.16 0.05 1.17 4.29 0.04 0.74 2.80 0.03 GMM 2.23 6.40 0.01 1.13 3.35 0.11 0.69 2.30 0.01

better mAW and mCCE scores than their Norm variants, which means the predicted intervals are more informative and better calibrated. On PEMSBay, the GMM variants also have better mAW, but Norm variants have better mCCE scores. This observation suggests that the underlying data distribution for PEMS-Bay is generally closer to a unimodal Gaussian. For the other two datasets, traffic dynamics are more complex and thus beyond the expressive power of unimodal Gaussians. Table 4 presents the deterministic performance metrics, which is a common practice in the literature (Li et al., 2018). The metrics include Mean Absolute Error (MAE), Mean Absolute Percentage Error (MAPE) and Root Mean Square Error (RMSE). When comparing across base models, STAEformer achieves the best deterministic performance on all three datasets, 18

Table 3: Relative CRPS improvement (%) of GMM over deterministic variants

Dataset

LGC STID GWNet MTGNN STAE Avg.

METR-LA 27.65 28.10 26.06 26.95 PEMS-Bay SimBarcaSpd 28.57 28.30

27.80 27.85 26.92

27.24 26.75 28.00

26.64 27.49 27.56 27.03 29.59 28.28

Average

27.52

27.33

27.93 27.60

27.43 27.78

which is consistent with the literature (Liu et al., 2023). In general, the GMM variants achieve better RMSE, and the natively deterministic models have better MAE and MAPE. This can be attributed to the fact that the deterministic models are trained with MAE as the loss function, and the GMM NLL loss in Equation 5 has an L2 form closer to the Mean Square Error (the square of RMSE). Since our experiments include no change in training pipelines, hyper-parameters, or model architectures, it is reasonable that the GMM variants perform better on RMSE while the deterministic ones are better on MAE and MAPE. 4.3. Performance Analysis by Prediction Horizons Figure 6 shows the CRPS scores of all methods at all prediction horizons, providing a more detailed view of the CRPS results in Table 2. A clear trend is that CRPS increases with prediction horizon for all methods, as predictions become more challenging. From the figures, a clear gap can be identified between the deterministic and probabilistic variants, which again demonstrates the drawbacks of false confidence and the necessity of probabilistic modeling in traffic forecasting. The performance hierarchy is mostly consistent across all prediction horizons, where a stronger backbone model (e.g. STAEformer) consistently outperforms a weaker one (e.g. STID), and the GMM variant consistently outperforms the Normal variant. Figure 7 shows the mAW scores of all probabilistic methods, which indicate the level of predicted uncertainty. As expected, the predicted uncertainty also generally increases with the prediction horizon, since longerhorizon predictions are more challenging. Notably, compared to the Normal variants, the GMM variants consistently have considerably lower mAW scores, indicating that the predictions are more informative with narrower confidence intervals. 19

Table 4: Deterministic performance metrics. Lower values are preferred for all, best result for a backbone in bold.

RMSE

MAPE

MAE

RMSE

MAPE

MAE

SimBarcaSpd

Det 3.11 8.69 6.20 1.65 3.92 3.60 1.05 24.06 1.76 Norm 3.35 9.42 6.13 1.78 4.33 3.69 1.12 25.91 1.71 GMM 3.24 9.08 6.01 1.73 4.30 3.68 1.10 25.93 1.70 Det 3.31 9.59 6.68 1.67 3.81 3.70 1.06 24.10 1.77 Norm 3.55 10.15 6.46 1.77 4.04 3.69 1.12 25.70 1.73 GMM 3.46 9.88 6.37 1.73 3.97 3.64 1.11 25.59 1.72

STAE MTGNN GWNet

LGC

RMSE

PEMS-Bay

STID

MAPE

METR-LA MAE

Method

Det 3.13 8.70 6.27 1.58 3.59 3.52 1.04 24.00 1.74 Norm 3.38 9.42 6.17 1.67 3.79 3.50 1.13 26.07 1.71 GMM 3.28 9.14 6.09 1.61 3.65 3.44 1.11 25.81 1.71 Det 3.12 8.64 6.37 1.57 3.53 3.52 1.00 22.91 1.69 Norm 3.40 9.36 6.18 1.66 3.77 3.46 1.08 25.14 1.65 GMM 3.28 9.08 6.10 1.62 3.69 3.45 1.06 24.63 1.65 Det 3.04 8.36 6.13 1.56 3.45 3.51 0.98 22.81 1.67 Norm 3.19 8.74 5.97 1.61 3.64 3.42 1.02 23.91 1.59 GMM 3.12 8.43 5.97 1.58 3.50 3.37 1.02 23.66 1.59

Figure 8 shows the calibration curves, i.e., the actual data coverage at different confidence levels. The difference between the curves and the diagonal line indicates the calibration quality, and is quantified by the mCCE metric in Equation 15. The calibration curve provides more details than the overall mCCE metric. In the ideal case, the actual data coverage should be perfectly aligned with the confidence levels. Any point above the ideal coverage line is under-confident, because the intervals cover more data than expected. On the contrary, any point below the ideal line is over-confident. Comparing across three datasets, complex backbone models generally have lower data coverage than simpler models. For example, STAEformer-based GMM and Normal variants have lower data coverage than corresponding STID-based variants on all three datasets. Meanwhile, GMM variants generally have lower data 20

LGC (Det) LGC (Norm) LGC (GMM) GWNet (Det) GWNet (Norm) GWNet (GMM) STID (Det) STID (Norm) STID (GMM) STAEformer (Det) STAEformer (Norm) STAEformer (GMM) MTGNN (Det) MTGNN (Norm) MTGNN (GMM)

CRPS

3.0

2.5

LGC (Det) LGC (Norm) LGC (GMM) GWNet (Det) GWNet (Norm) GWNet (GMM) STID (Det) STID (Norm) STID (GMM) STAEformer (Det) STAEformer (Norm) STAEformer (GMM) MTGNN (Det) MTGNN (Norm) MTGNN (GMM)

2.0 1.8 1.6 CRPS

3.5

1.4 1.2 1.0

2.0

1.1 1.0

0.7

0.6 2

3

4

5

6 7 8 Prediction Horizon

9

10

11

0.9 0.8

0.8 1

LGC (Det) LGC (Norm) LGC (GMM) GWNet (Det) GWNet (Norm) GWNet (GMM) STID (Det) STID (Norm) STID (GMM) STAEformer (Det) STAEformer (Norm) STAEformer (GMM) MTGNN (Det) MTGNN (Norm) MTGNN (GMM)

1.2

CRPS

4.0

1

12

2

(a) METR-LA

3

4

5

6 7 8 Prediction Horizon

9

10

11

0.6

12

1

(b) PEMS-Bay

2

3

4

5 6 7 Prediction Horizon

8

9

10

(c) SimBarcaSpd

Figure 6: CRPS score of all methods at all prediction horizons on three datasets.

10 9

6 Average Width (AW) of CI

Average Width (AW) of CI

11

8 7

3.8

LGC (Norm) LGC (GMM) GWNet (Norm) GWNet (GMM) STID (Norm) STID (GMM) STAEformer (Norm) STAEformer (GMM) MTGNN (Norm) MTGNN (GMM)

5

3.4

4

3

1

2

3

4

5

6 7 8 Prediction Horizon

9

(a) METR-LA

10

11

12

3.2 3.0 2.8 2.6 2.4

6

2

LGC (Norm) LGC (GMM) GWNet (Norm) GWNet (GMM) STID (Norm) STID (GMM) STAEformer (Norm) STAEformer (GMM) MTGNN (Norm) MTGNN (GMM)

3.6

Average Width (AW) of CI

LGC (Norm) LGC (GMM) GWNet (Norm) GWNet (GMM) STID (Norm) STID (GMM) STAEformer (Norm) STAEformer (GMM) MTGNN (Norm) MTGNN (GMM)

12

2.2

1

2

3

4

5

6 7 8 Prediction Horizon

9

(b) PEMS-Bay

10

11

12

1

2

3

4

5 6 7 Prediction Horizon

8

9

10

(c) SimBarcaSpd

Figure 7: mAW score of all probabilistic methods at all prediction horizons on three datasets.

coverage than Normal variants at the same confidence levels, indicating that the GMM confidence intervals are more selective. On METR-LA and SimBarcaSpd, the GMM variants of all models are notably closer to the ideal calibration than Normal variants, which is also indicated by the mCCE metric in Table 2. However, on PEMS-Bay, the Normal variants generally have better calibration than GMM variants. This difference suggests that the underlying traffic in PEMS-Bay is closer to a single-mode Gaussian distribution, while METR-LA and SimBarcaSpd have more complex multi-modal distributions that can be better captured by GMMs. Such an observation is consistent with the deterministic performance in Table 4, where the percentage error (MAPE) on PEMS-Bay is around 4%, meaning that the traffic speed in this region stays within a predictable range. In contrast, all methods have significantly higher MAPE on METR-LA (around 9%) and SimBarcaSpd (around 24%), indicating more volatile traffic conditions. 21

0.9

0.9

0.9

Ideal Coverage LGC (Norm) LGC (GMM) GWNet (Norm) GWNet (GMM) STID (Norm) STID (GMM) STAEformer (Norm) STAEformer (GMM) MTGNN (Norm) MTGNN (GMM)

0.7

0.6

0.5 0.5

0.6

0.7 Confidence

0.8

0.7

Ideal Coverage LGC (Norm) LGC (GMM) GWNet (Norm) GWNet (GMM) STID (Norm) STID (GMM) STAEformer (Norm) STAEformer (GMM) MTGNN (Norm) MTGNN (GMM)

0.6 0.5 0.4

0.9

0.5

0.6

(a) METR-LA

0.7 Confidence

0.8

Data Coverage

0.8

0.8

Data Coverage

Data Coverage

0.8

Ideal Coverage LGC (Norm) LGC (GMM) GWNet (Norm) GWNet (GMM) STID (Norm) STID (GMM) STAEformer (Norm) STAEformer (GMM) MTGNN (Norm) MTGNN (GMM)

0.7

0.6

0.5

0.9

0.5

(b) PEMS-Bay

0.6

0.7 Confidence

0.8

0.9

(c) SimBarcaSpd

Figure 8: Calibration curve of all probabilistic methods at different confidence interval levels on three datasets.

8

6

5

4 2

4

6 8 Prediction Horizon

(a) METR-LA

10

12

4.0

LGC (Det) LGC (Norm) LGC (GMM) GWNet (Det) GWNet (Norm) GWNet (GMM) STID (Det) STID (Norm) STID (GMM) STAEformer (Det) STAEformer (Norm) STAEformer (GMM) MTGNN (Det) MTGNN (Norm) MTGNN (GMM)

3.5 3.0 2.5 2.0 1.5 2

4

6 8 Prediction Horizon

(b) PEMS-Bay

10

12

1.9 1.8 rmse

LGC (Det) LGC (Norm) LGC (GMM) GWNet (Det) GWNet (Norm) GWNet (GMM) STID (Det) STID (Norm) STID (GMM) STAEformer (Det) STAEformer (Norm) STAEformer (GMM) MTGNN (Det) MTGNN (Norm) MTGNN (GMM)

rmse

7

rmse

2.0

4.5

1.7 1.6

LGC (Det) LGC (Norm) LGC (GMM) GWNet (Det) GWNet (Norm) GWNet (GMM) STID (Det) STID (Norm) STID (GMM) STAEformer (Det) STAEformer (Norm) STAEformer (GMM) MTGNN (Det) MTGNN (Norm) MTGNN (GMM)

1.5 2

4

6 Prediction Horizon

8

10

(c) SimBarcaSpd

Figure 9: RMSE of all methods at all prediction horizons on three datasets.

Figure 9 shows the RMSE of all methods at all prediction horizons, which increases with the prediction horizon as expected. Still, the shape of the curves in SimBarcaSpd differs from the other two datasets, where RMSE increases more rapidly at the last few horizons. Unlike the highway traffic in METR-LA and PEMS-Bay, SimBarcaSpd comes from an urban scenario with dense connections between road segments, and predicting longer horizons therefore becomes increasingly challenging due to complex interactions. 4.4. Prediction Examples To further illustrate the challenge of modeling urban traffic, Figure 10 shows the test-set data distribution for three representative road segments in SimBarcaSpd. From these visualizations, it is clear that the traffic speed for urban roads can have multiple modes (e.g., congested vs free-flow), while the middle range is less likely. Figure 11 presents the 30-minute-ahead prediction examples of the LGC GMM at these segments; data from a high-traffic22

12

10

10

8 6 4

8 6 4

2

2

0

0 0

20

40 60 80 Time from Simulation Start (min)

(a) Segment 113

100

10 Segment Speed (m/s)

14

12 Segment Speed (m/s)

Segment Speed (m/s)

14

8 6 4 2 0

0

20

40 60 80 Time from Simulation Start (min)

(b) Segment 345

100

0

20

40 60 80 Time from Simulation Start (min)

100

(c) Segment 1068

Figure 10: Test set data distribution in SimBarcaSpd. The data from 26 test set sessions are plotted (colored by session) according to their timestamp in a day. The scatter plot clearly shows that the traffic speed at these segments has multiple modes, e.g., a congested mode and a free-flow mode.

demand scenario are selected for visualization. In both segments 113 and 345, the deterministic prediction is a stable-looking curve around a middle value, yet road speeds in reality flip between low and high values. Under such scenarios, deterministic predictions are very misleading, but this is a fundamental limitation coming from the training paradigm, where the model is forced to predict a single value to minimize the absolute error. On the contrary, our GMM adaptation allows the model to produce multi-modal distributions, and it is clear that for both segments, the GMM predictions have two distinct modes that can cover the two possibilities in reality. For segment 345, the model can predict a high-speed mode around 11 m/s and a low-speed one around 2 m/s, both aligned with the scatter points in the data distribution as shown in Figure 10. Another advantage of multi-modal predictions is faster adaptation to changes. For segment 1068, traffic congestion builds up from 75 min, and the deterministic prediction reacts slowly, while the mode for the congested regime appears in the GMM prediction as early as 78 min. This ability can be crucial in traffic management, where earlier warnings of congestion allow more effective mitigation measures. Figure 12 shows how the prediction time horizon affects the predicted distributions. For the same simulation time (e.g., 60 min to 90 min), the onestep-ahead prediction has a sharper single-mode distribution, since predicting the immediate next step is relatively easier. As the horizon extends to five and ten steps, the predicted distribution becomes more uncertain for the same time span, with wider spread and multiple modes. This behavior shows the GMM model’s ability to capture the increasing uncertainty with longer future horizons. 23

14

30min Pred Ground Truth

12

6

10 Speed (m/s)

Speed (m/s)

8

8 6

8 6

4

4

2

2

2

0

0

(a) Segment 113

60 63 66 69 72 75 78 81 84 87 90 93 96 99 102 105 108 111 114 117

4

Time from Simulation Start (min)

30min Pred Ground Truth

12

10

60 63 66 69 72 75 78 81 84 87 90 93 96 99 102 105 108 111 114 117

Speed (m/s)

10

14

30min Pred Ground Truth

12

Time from Simulation Start (min)

(b) Segment 345

0

60 63 66 69 72 75 78 81 84 87 90 93 96 99 102 105 108 111 114 117

14

Time from Simulation Start (min)

(c) Segment 1068

Figure 11: The 30-minute-ahead predictions (10 steps) of LGC (GMM) for a high-trafficdemand scenario in SimBarcaSpd. The density ridges are colored for better visualization.

Figure 12: One, five, and ten step predictions of LGC (GMM) for road segment 806 in SimBarcaSpd. The density ridges are colored for better visualization.

4.5. Robustness to Data Quality To investigate the impact of data quality on probabilistic traffic forecasting, we conduct an additional experiment using the SimBarca dataset (Xiong et al., 2025), featuring drone-measured high-frequency traffic speed data collected every 5 seconds. Concretely, SimBarca and our derived SimBarcaSpd dataset have the same input and output window (30 min), but SimBarca has finer temporal resolution (5 seconds), providing more detailed traffic dynamics than SimBarcaSpd (3 minutes). We adapt the LGC model to SimBarca by adding a temporal convolution layer before the LSTM to downsample the higher-resolution input, while keeping the rest of the architecture unchanged. We study the impact of data quality from two aspects: • Partial coverage: when sensors cover only a random 10% of the road segments (Xiong et al., 2025). • Temporal detail: when input data are provided every 5 seconds (more details) or 3 minutes (marked as Low Res). 24

Table 5: Average CRPS with different input data quality. Numbers in parentheses indicate the percentage increase compared to the ideal case, lower values are preferred.

Model Ideal

10% coverage

Low resolution

GMM 0.7505 0.7770 (+3.5%) 0.7547 (+0.6%) Det 1.0104 1.0858 (+7.5%) 1.0543 (+4.3%) Figure 13 shows the CRPS scores of the LGC model on SimBarcaSpd (marked as Low Res) and the adapted LGC model on SimBarca under different data quality settings. Table 5 also summarizes the average CRPS across all prediction horizons for quantitative comparison. The percentage in parentheses indicates the increased CRPS compared to the ideal case. The ideal case is when the model has access to high-resolution data with full coverage, and both models perform best compared with degraded dataquality settings. When only 10% of the road segments are covered, the CRPS for the probabilistic GMM method increases by 3.5%, while the deterministic Det method has a larger CRPS increase of 7.5%. This observation is also consistent with the larger gap in Figure 13 between the ideal and 10% coverage curves for Det compared to GMM. The impact of lower temporal resolution is less significant than coverage, where the CRPS for GMM increases by only 0.6% and Det increases by 4.3%. These comparisons suggest that multi-modal probabilistic modeling can provide more accurate and informative predictions while also being more robust to data-quality issues such as sensor coverage and data frequency. 5. Conclusions In summary, this paper presents a universal approach to transform existing spatio-temporal forecasting models into probabilistic ones by replacing only the final output layer with a novel GMM layer. Our adaptation method is simple, effective, and can be applied to a wide range of backbone architectures with minimal modifications to model structure and training procedure. A comprehensive evaluation methodology has been introduced to assess the quality of predictive distributions, through which we show that probabilistic modeling consistently outperforms deterministic approaches across multiple datasets and architectures. Compared to deterministic predictions, which are often vulnerable to false-confidence issues, our multi-modal probabilistic 25

100%-GMM 100%-Det 100%-GMM (Low Res) 100%-Det (Low Res) 10%-GMM 10%-Det

1.2 1.1

CRPS

1.0 0.9 0.8 0.7 1

2

3

4

5 6 7 Prediction Horizon

8

9

10

Figure 13: CRPS scores of LGC variants on SimBarca under different data quality settings.

predictions can better capture the complex and multi-modal nature of traffic data, providing more accurate, more informative, and safer predictions for real-world applications. Additional analyses have demonstrated that a probabilistic model equipped with the GMM prediction layer can be more robust to data quality issues such as coverage (missing rate) and temporal resolution. We hope our work can raise awareness of probabilistic modeling in the community and provide a novel perspective for evaluating spatio-temporal forecasting models. Future research could explore the potential of multimodal distributions with more advanced base architectures and probabilistic modeling techniques. CRediT Author Contributions Weijiang Xiong: Conceptualization, Methodology, Data Curation, Software, Writing - Original Draft. Robert Fonod: Conceptualization, Methodology, Writing - Review & Editing, Supervision. Nikolas Geroliminis: Conceptualization, Methodology, Writing - Review & Editing, Supervision, Funding Acquisition. Acknowledgements This work is supported by the Swiss National Science Foundation (SNSF) under NCCR Automation, grant agreement number 51NF40_180545. 26

References An, Y., Li, Z., Li, X., Liu, W., Yang, X., Sun, H., Chen, M., Zheng, Y., Gong, Y., 2025. Spatio-temporal multivariate probabilistic modeling for traffic prediction. IEEE Transactions on Knowledge and Data Engineering 37, 2986–3000. Bessa, R.J., Trindade, A., Miranda, V., 2014. Spatial-temporal solar power forecasting for smart grids. IEEE Transactions on Industrial Informatics 11, 232–241. Candille, G., Talagrand, O., 2005. Evaluation of probabilistic prediction systems for a scalar variable. Quarterly Journal of the Royal Meteorological Society: A journal of the atmospheric sciences, applied meteorology and physical oceanography 131, 2131–2150. Chen, X., Chen, R., 2019. A review on traffic prediction methods for intelligent transportation system in smart cities, in: 2019 12th International Congress on Image and Signal Processing, BioMedical Engineering and Informatics (CISP-BMEI), IEEE. pp. 1–5. Coric, V., Wang, Z., Vucetic, S., 2011. Traffic speed forecasting by mixture of experts, in: 2011 14th International IEEE Conference on Intelligent Transportation Systems, IEEE. Ermagun, A., Levinson, D., 2018. Spatiotemporal traffic forecasting: review and proposed directions. Transport Reviews 38, 786–814. Ferreira, G.O., Ravazzi, C., Dabbene, F., Calafiore, G.C., Fiore, M., 2023. Forecasting network traffic: A survey and tutorial with open-source comparative evaluation. IEEE Access 11, 6018–6044. Fu, J., Zhou, W., Chen, Z., 2021. Bayesian spatio-temporal graph convolutional network for traffic forecasting, in: Proceedings of the AAAI Conference on Artificial Intelligence, AAAI. Gal, Y., Ghahramani, Z., 2016. Dropout as a bayesian approximation: Representing model uncertainty in deep learning, in: international conference on machine learning, PMLR. pp. 1050–1059. Gelman, A., Carlin, J.B., Stern, H.S., Rubin, D.B., 1995. Bayesian data analysis. Chapman and Hall/CRC. 27

Ghosh, B., Basu, B., O’Mahony, M., 2010. Random process model for urban traffic flow using a wavelet-bayesian hierarchical technique. ComputerAided Civil and Infrastructure Engineering 25, 613–624. Gneiting, T., Raftery, A.E., 2007. Strictly proper scoring rules, prediction, and estimation. Journal of the American statistical Association 102, 359– 378. Guo, C., Pleiss, G., Sun, Y., Weinberger, K.Q., 2017. On calibration of modern neural networks, in: International conference on machine learning, PMLR. pp. 1321–1330. Guo, S., Lin, Y., Feng, N., Song, C., Wan, H., 2019. Attention based spatialtemporal graph convolutional networks for traffic flow forecasting, in: Proceedings of the AAAI conference on artificial intelligence, pp. 922–929. Hochreiter, S., Schmidhuber, J., 1997. Long short-term memory. Neural Computation 9, 1735–1780. Jiang, W., Luo, J., 2022. Graph neural network for traffic forecasting: A survey. Expert systems with applications 207, 117921. Kendall, A., Gal, Y., 2017. What uncertainties do we need in bayesian deep learning for computer vision? Advances in neural information processing systems 30. Kipf, T.N., Welling, M., 2016. Semi-supervised classification with graph convolutional networks. arXiv preprint arXiv:1609.02907 . Kwak, S., Geroliminis, N., Frossard, P., 2021. Traffic signal prediction on transportation networks using spatio-temporal correlations on graphs. IEEE Transactions on Signal and Information Processing over Networks 7, 648–659. doi:10.1109/TSIPN.2021.3118489. Li, R., Zhang, F., Li, T., Zhang, N., Zhang, T., 2022. Dmgan: Dynamic multi-hop graph attention network for traffic forecasting. IEEE Transactions on Knowledge and Data Engineering 35, 9088–9101. Li, Y., Yu, R., Shahabi, C., Liu, Y., 2018. Diffusion convolutional recurrent neural network: Data-driven traffic forecasting, in: International Conference on Learning Representations (ICLR ’18). 28

Lim, B., Zohren, S., 2021. Time-series forecasting with deep learning: a survey. Philosophical transactions of the royal society a: mathematical, physical and engineering sciences 379. Liu, H., Dong, Z., Jiang, R., Deng, J., Deng, J., Chen, Q., Song, X., 2023. Spatio-temporal adaptive embedding makes vanilla transformer sota for traffic forecasting, in: Proceedings of the 32nd ACM international conference on information and knowledge management, pp. 4125–4129. Liu, J., Wang, Y., Zhu, J., Bai, W., Zhang, H., Zuo, L., Zhou, T., Li, K., 2025. A multilayer spatiotemporal correlation-aware graph attention network for traffic flow prediction. IEEE Transactions on Neural Networks and Learning Systems . Loshchilov, I., Hutter, F., 2019. Decoupled weight decay regularization, in: International Conference on Learning Representations. URL: https:// openreview.net/forum?id=Bkg6RiCqY7. Maclean, I.M., 2020. Predicting future climate at high spatial and temporal resolution. Global change biology 26, 1003–1011. Mauro, R., et al., 2015. Traffic and random processes. Trento: Springer International Publishing Switzerland. DOI: https://doi. org/10.1007/9783-319-09324-6 . Qian, W., Nielsen, T.D., Zhao, Y., Larsen, K.G., Yu, J.J., 2024. Uncertaintyaware temporal graph convolutional network for traffic speed forecasting. IEEE Transactions on Intelligent Transportation Systems 25, 8578–8591. Qian, W., Zhao, Y., Zhang, D., Chen, B., Zheng, K., Zhou, X., 2023. Towards a unified understanding of uncertainty quantification in traffic flow forecasting. IEEE Transactions on Knowledge and Data Engineering 36, 2239–2256. Ramezani, M., Geroliminis, N., 2012. On the estimation of arterial route travel time distribution with markov chains. Transportation Research Part B: Methodological 46, 1576–1590. Ravish, R., Swamy, S.R., 2021. Intelligent traffic management: A review of challenges, solutions, and future perspectives. Transport and Telecommunication 22, 163–182. 29

Sengupta, A., Mondal, S., Das, A., Guler, S.I., 2024. A bayesian approach to quantifying uncertainties and improving generalizability in traffic prediction models. Transportation Research Part C: Emerging Technologies 162, 104585. Shao, Z., Zhang, Z., Wang, F., Wei, W., Xu, Y., 2022. Spatial-temporal identity: A simple yet effective baseline for multivariate time series forecasting, in: Proceedings of the 31st ACM International Conference on Information & Knowledge Management, pp. 4454–4458. Tan, Q., Ye, M., Ma, A.J., Yang, B., Yip, T.C.F., Wong, G.L.H., Yuen, P.C., 2020. Explainable uncertainty-aware convolutional recurrent neural network for irregular medical time series. IEEE Transactions on Neural Networks and Learning Systems 32, 4665–4679. Tang, B., Matteson, D.S., 2021. Probabilistic transformer for time series analysis. Advances in neural information processing systems 34, 23592– 23608. Tedjopurnomo, D.A., Bao, Z., Zheng, B., Choudhury, F.M., Qin, A.K., 2022. A survey on modern deep neural network for traffic prediction: Trends, methods and challenges. IEEE Transactions on Knowledge and Data Engineering 34, 1544–1561. doi:10.1109/TKDE.2020.3001195. Tsitsokas, D., Kouvelas, A., Geroliminis, N., 2023. Two-layer adaptive signal control framework for large-scale dynamically-congested networks: Combining efficient max pressure with perimeter control. Transportation Research Part C: Emerging Technologies 152, 104128. Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A.N., Kaiser, Ł., Polosukhin, I., 2017. Attention is all you need. Advances in Neural Information Processing Systems 30. Vlahogianni, E.I., Karlaftis, M.G., Golias, J.C., 2014. Short-term traffic forecasting: Where we are and where we’re going. Transportation Research Part C: Emerging Technologies 43, 3–19. Wang, Y., Yang, X., Liang, H., Liu, Y., 2018. A review of the self-adaptive traffic signal control system based on future traffic environment. Journal of Advanced Transportation 2018, 1096123. 30

Wanke, C., Mulgund, S., Greenbaum, D., Song, L., 2004. Modeling traffic prediction uncertainty for traffic management decision support, in: AIAA Guidance, Navigation, and Control Conference and Exhibit, p. 5230. Wu, Y., James, J., 2021. A bayesian learning network for traffic speed forecasting with uncertainty quantification, in: 2021 International Joint Conference on Neural Networks (IJCNN), IEEE. pp. 1–7. Wu, Y., Ye, Y., Zeb, A., Yu, J.J., Wang, Z., 2023. Adaptive modeling of uncertainties for traffic forecasting. arXiv preprint arXiv:2303.09273 . Wu, Z., Pan, S., Long, G., Jiang, J., Chang, X., Zhang, C., 2020. Connecting the dots: Multivariate time series forecasting with graph neural networks, in: Proceedings of the 26th ACM SIGKDD international conference on knowledge discovery & data mining, pp. 753–763. Wu, Z., Pan, S., Long, G., Jiang, J., Zhang, C., 2019. Graph wavenet for deep spatial-temporal graph modeling. arXiv preprint arXiv:1906.00121 . Xiao, X., Xiang, X., Yang, X., Jin, Z., Xu, J., Wang, S., Mao, G., Shao, W., 2026. Heuristic knowledge-driven spatio-temporal forecasting via multigraph. IEEE Transactions on Neural Networks and Learning Systems . Xiong, W., Fonod, R., Alahi, A., Geroliminis, N., 2025. Multi-source urban traffic flow forecasting with drone and loop detector data. IEEE Transactions on Intelligent Transportation Systems 26, 18637–18652. doi:10. 1109/TITS.2025.3605014. Yu, B., Yin, H., Zhu, Z., 2018. Spatio-temporal graph convolutional networks: A deep learning framework for traffic forecasting, in: Proceedings of the 27th International Joint Conference on Artificial Intelligence, AAAI Press. pp. 3634–3640. Yuan, H., Li, G., 2021. A survey of traffic prediction: from spatio-temporal data to intelligent transportation. Data Science and Engineering 6, 63–85. Zheng, C., Fan, X., Wang, C., Qi, J., 2020. GMAN: A graph multi-attention network for traffic prediction, in: Proceedings of the AAAI conference on artificial intelligence, pp. 1234–1241.

31

Zheng, Z., Zhang, Z., 2023. A stochastic recurrent encoder decoder network for multistep probabilistic wind power predictions. IEEE Transactions on Neural Networks and Learning Systems 35, 9565–9578. Zhou, Z., Gu, Z., Liu, P., Yu, W., Liu, Z., 2025. Leveraging semi-supervised learning and meta-learning for re-identification in few-shot spatiotemporal anomaly detection. IEEE Transactions on Neural Networks and Learning Systems .

32