Data-Driven Estimation of End-to-End Delay Probability Density Function for Time-Sensitive WiFi Networks

1. Introduction

Advancements in communication, computation, and artificial intelligence continuously promote the development of time-sensitive applications [1], such as industrial automation, intelligent transportation, telemedicine, data center networking, etc. Time-sensitive applications demand that the End-to-End (E2E) transmission delay remains within the millisecond range and that the delay jitter stays within the microsecond range, both with a high probability [2,3,4]. As depicted in Figure 1, for Ethernet, Time-Sensitive Networking (TSN) technology employs traffic shaping mechanisms such as IEEE 802.1Qav and IEEE 802.1Qbv to allocate deterministic transmission slots for high-priority traffic. This ensures the availability of transmission bandwidth and precise timing for critical operations, thereby enhancing the reliability and delay determinism of key services [5]. However, in wireless networks, factors such as dynamic node behavior, energy limitations, and perpetually evolving environments pose significant challenges to the implementation of TSN technologies. 5G ultra-reliable low-latency communication (URLLC) technology has achieved remarkable progress in reducing E2E delay. However, it does not adequately account for queueing delay, which may still pose challenges in certain real-time applications. As a result, there remains a gap in meeting the stringent delay requirements of time-sensitive applications. Consequently, ensuring a deterministic E2E delay remains a critical objective in the field of wireless networks. To achieve deterministic control of E2E delay in wireless networks, accurate delay estimation is essential, as relying solely on average delay is insufficient. Therefore, the delay probability density function (PDF) must be estimated so that resource and queue management strategies can be adjusted accordingly. This ensures the timely transmission of critical data even as network conditions change.

In the future, WiFi 7 is expected to incorporate TSN capabilities, thereby enabling ultra-reliability and low-latency in unlicensed frequency bands [2]. In contrast to mobile networks, the WiFi environment is predominantly indoor-oriented. Consequently, delay fluctuations in WiFi networks are less pronounced compared to those in mobile networks. Therefore, the estimation of the E2E delay PDF in WiFi networks is not as complex as that in mobile networks, yet similarities exist. In this paper, the primary focus is on the estimation of the E2E delay PDF in WiFi networks. It is expected that the associated findings can be extended and applied to mobile networks. At present, the estimation methods for network delay PDFs can be categorized into two types: theoretical analysis approaches and data-driven methods. Theoretical analysis methods generally depend on a set of assumptions and utilize mathematical tools to solve and analyze the system behavior. However, these assumptions may not always align with practical scenarios. In recent years, the integration of machine learning into network delay estimation has garnered significant attention [6]. The development in this field is primarily driven by data-driven approaches, where machine learning techniques based on neural networks are utilized to identify and learn the relationships between delay and other variables by extracting insights from real-world data.

As a representative data-driven approach, Mixture Density Networks (MDNs) have been extensively applied and validated for delay PDF estimation in complex systems. MDNs integrate neural networks with probabilistic models, such as the Gaussian Mixture Model (GMM), to estimate the parameters of the mixture model via a fully connected neural network, thereby generating a conditional PDF. Recently, Mostafavi et al. [7] employed MDN in conjunction with the Generalized Pareto Distribution (GPD) to estimate the E2E delay PDF of several real networks, including Commercial Off-The-Shelf (COTS) 5G, Open Air Interface (OAI) 5G, and Mango COMM IEEE 802.11 (WiFi). MDN maps the transmission conditions or traffic characteristics to the parameters of the PDF and generates the PDF. However, for WiFi networks, there is a noticeable gap between the estimation results and actual measurements, particularly in the tail of the PDF. To address this issue, this paper proposes a two-stage and three-divided estimation method, referred to as TSTD, for estimating the E2E delay PDF in WiFi networks. The main contributions are summarized as follows.

• A three-divided PDF is constructed, consisting of a main part and a two-segment tail part. The main part is modeled using GMM, while the two segments of the tail part are modeled using distinct GPDs. Additionally, the coupling relationships between the different segments in the main part and the tail part are taken into account.
• For the parameters of the PDF, the thresholds of the GPDs are determined based on the variation trend of the tail probability in E2E delay measurements, and the remaining parameters are obtained through training via MDN.
• Compared with existing methods, the E2E delay PDF obtained using the proposed method exhibits a smaller discrepancy from actual measurements. In particular, the mean absolute errors and average fluctuation amplitudes of tail probabilities at certain delay values are reduced by at least one order of magnitude. Moreover, the multiple-segmentation feature of the proposed method enhances its robustness in situations where measurement data are affected by low levels of Gaussian noise.

The remainder of this paper is structured as follows. Section 2 reviews the related work. Section 3 provides a concise introduction to the research problem. Section 4 elaborates on the proposed two-stage and three-divided estimation method. Section 5 details the experimental procedures and performance analysis. Finally, Section 6 presents the conclusion and discusses future work.

2. Related Work

The estimation methods for the E2E delay probability distribution in multiple application scenarios have garnered significant attention. Cao et al. [8] modeled the Mobile Edge Computing (MEC) network as a two-stage tandem queueing system. Subsequently, they proposed an estimation method for the probability distribution of the E2E delay using the matrix-geometric method. Mei et al. [9] modeled and analyzed the delay bound of a multi-cluster MEC network using the stochastic network calculus (SNC) approach. Cui et al. [10] proposed a method for calculating the E2E delay violation probability of target traffic in the industrial Internet of Things (IoT), also based on SNC. Coll-Perales et al. [11] proposed a 5G E2E delay model for Vehicle-to-Network (V2N) and Vehicle-to-Network-to-Vehicle (V2NV2N) communications. They quantified and analyzed the 5G E2E delay performance for V2N and V2NV2N communications under various 5G network deployments and configurations.

The aforementioned methods are all theoretical analysis. Theoretical analysis methods generally rely on a series of assumptions and employ mathematical tools to solve and analyze the system’s behavior. However, these assumptions may not always align with practical scenarios. In view of this, data-driven methods have been extensively studied for estimating network delay probability distributions.

Fadhil et al. [12] proposed a method for modeling the E2E delay of 5G networks using a GMM. Based on E2E delay data, an Expectation-Maximization (EM) algorithm was employed to estimate the GMM parameters. The results indicate that as the number of data samples and GMM components increases, higher accuracy can be achieved. However, the computation time also increases rapidly with the increase in the number of samples and components. Specifically, the increasing trend is approximated as an exponential behavior as the number of GMM components increases [12]. In [13], Chen et al. employed the deep learning approach to generate the probability density function (PDF) of practical data. The cumulative distribution functions (CDFs) of common probability distributions were utilized as activation functions in the hidden layers of the proposed deep learning model to learn actual cumulative probabilities. Furthermore, the differential equation derived from the trained deep learning model can be used to estimate the PDF. Experimental results demonstrated that both the CDF and PDF can be accurately estimated by the proposed method, as assessed by the mean absolute percentage error. In [14], the MDN combined with GMM was utilized to estimate the delay probability distribution of single-stage queueing systems. Subsequently, this work was extended by Raeis et al. [15] to more complex systems, namely service function chains. In [14,15], the estimation methods have been proven effective in fitting the main part of the delay probability distribution using machine learning techniques. However, the fitting results for the tail part of the distribution are less satisfactory. This is because the output of GMM-based models typically exhibits exponentially decaying tail probabilities, which can lead to significant errors when dealing with low-probability events, especially in the case of heavy-tailed distributions.

It is shown in some literature that better fitting results can be obtained when the main and tail parts of the delay probability distribution are fitted using different functions respectively. Yasuda et al. [16] utilized the bistate model to simulate wireless networks involving connection and disconnection processes. For the E2E delay probability distribution, the main part, which arises during the connection process, was fitted using the shifted Gamma distribution. In contrast, the tail part, resulting from the accumulation of probe packets during the disconnection phase, was approximated by the exponential distribution. This approach exhibited superior performance compared to traditional methods in terms of both the negative log predictive density (NLPD) [17] and the continuous ranked probability score (CRPS) [18]. In addition, the GPD is a classical asymptotically motivated model for the unknown excess distribution above high thresholds [19,20]. Consequently, it has been widely applied to fit the tail part of the delay PDF [7,21]. In [21], Mostafavi et al. proposed an extreme value mixture model based on the mixture density network (MDN). This was achieved by integrating the GPD tail model with the GMM. Specifically, the GMM was used to capture the main part of the delay PDF, while the GPD was employed to characterize the tail of the delay PDF. Numerical experiments conducted on a three-stage tandem queueing system demonstrated that the proposed method outperforms existing state-of-the-art GMM-based estimation techniques. Furthermore, Mostafavi et al. [7] integrated the MDN with the GPD to estimate the PDF of E2E delay in three wireless network scenarios, namely commercial off-the-shelf (COTS) 5G, Open Air Interface (OAI) 5G, and Mango COMM IEEE 802.11 (WiFi). The GMM was employed to fit the main part of the delay probability distribution, while the tail part was modeled using the GPD. For the 5G scenarios, this approach exhibited robust performance through noise regularization when the tail profile was nonsmooth. However, in the WiFi network scenario, the estimation results showed a significant gap and fluctuation compared to actual measurements.

In [7], the threshold parameter of the GPD is obtained through training. It is worth noting that the threshold parameter of the GPD can also be determined using other methods, such as empirical approaches and graphical techniques. The empirical method depends on an analyst’s statistical knowledge and expertise to select an appropriate threshold, which may introduce bias [22,23]. Graphical diagnostic methods [24] have been widely used for threshold estimation, utilizing specific plots of the measured data to observe trends and assist in threshold selection. Cyrille et al. [25] utilized graphical diagnostic tools to determine the threshold range. Based on this, they established specific thresholds and optimized the parameters of the GPD using the Kolmogorov-Smirnov (KS) goodness-of-fit test. They further validated that the estimated parameters exhibited improved accuracy with larger sample sizes through the use of actual hydrological data. Zhao et al. [26] pointed out that the process of GPD parameter estimation is independent of threshold selection. They compared the performance of three threshold selection procedures to determine an appropriate threshold for asymptotic fitting of the GPD above this threshold. The results demonstrated that this estimator achieved satisfactory performance in environmental data analysis, particularly in fitting the tail probability distribution.

Drawing on the strengths of the existing works mentioned above, this paper proposes an improved estimation method to enhance the accuracy of fitting the E2E delay PDF in WiFi networks. Firstly, the threshold parameter of the GPD is determined based on the variation trend of E2E delay measurements. Subsequently, the remaining parameters of the GPD are estimated via the MDN. Moreover, the tail part of the delay PDF is divided into two segments, each of which is fitted using a distinct GPD.

3. Research Problem Statement

The conditional PDF for the E2E delay Y of packets with length X in WiFi networks is estimated based on the measured delay data, as illustrated in Figure 2.

Let $\left\{(x_n,y_n):n=0,1,\dots ,N\right\}$ denote the delay dataset measured over a period of time, where $(x_n,y_n)$ is the record corresponding to the n-th packet, $x_n$ denotes the packet length, $x_n\in S=\left\{s_k:k=0,1,\dots ,K\right\}$, $s_k>0$, and $y_n(y_n>0)$ denotes the E2E (uplink or downlink) delay. For packets with length $X=x\in S$, the conditional PDF $h\left(y|\widehat{\theta }\left(x\right)\right)$ will be estimated using some samples from the delay dataset, where $\widehat{\theta }\left(x\right)$ represents the estimated parameter vector. In additional, the CDF is denoted as $H\left(y|\widehat{\theta }\left(x\right)\right)=P\left\{Y\le y|\widehat{\theta }\left(x\right)\right\}$. The tail probability is denoted as $\overline{H}\left(y|\widehat{\theta }\left(x\right)\right)=1-H\left(y|\widehat{\theta }\left(x\right)\right)$.

4. Two-Stage and Three-Divided Estimation Method

The proposed method, TSTD, consists of two steps: namely, PDF construction and parameter estimation, as described below.

$$h\left(y\mid \widehat{\theta }\left(x\right)\right)=\left\{\begin{array}{cc}f\left(y\mid {\widehat{\theta }}_0\left(x\right)\right),\hfill & y\le u_1\left(x\right),\hfill \\ {\alpha }_1\left(x\right)\left[1-F\left(u_1\left(x\right)\mid {\widehat{\theta }}_0\left(x\right)\right)\right]g_1\left(y\mid {\widehat{\theta }}_1\left(x\right)\right)+\left(1-{\alpha }_1\left(x\right)\right)f\left(y\mid {\widehat{\theta }}_0\left(x\right)\right),\hfill & u_1\left(x\right)<y\le u_2\left(x\right),\hfill \\ {\alpha }_1\left(x\right){\alpha }_2\left(x\right)\left[1-F\left(u_1\left(x\right)\mid {\widehat{\theta }}_0\left(x\right)\right)\right]\left[1-G_1\left(u_2\left(x\right)\mid {\widehat{\theta }}_1\left(x\right)\right)\right]g_2\left(y\mid {\widehat{\theta }}_2\left(x\right)\right)\hfill \\ +{\alpha }_1\left(x\right)\left(1-{\alpha }_2\left(x\right)\right)\left[1-F\left(u_1\left(x\right)\mid {\widehat{\theta }}_0\left(x\right)\right)\right]g_1\left(y\mid {\widehat{\theta }}_1\left(x\right)\right)\hfill \\ +\left(1-{\alpha }_1\left(x\right)\right)f\left(y\mid {\widehat{\theta }}_0\left(x\right)\right),\hfill & y>u_2\left(x\right).\hfill \end{array}\right.$$

(1)

4.1. PDF construction

The PDF $h\left(y|\widehat{\theta }\left(x\right)\right)$ is expressed in the form of (1), where

$$\widehat{\theta }\left(x\right)=\left({\alpha }_1\left(x\right),{\alpha }_2\left(x\right),{\widehat{\theta }}_0\left(x\right),{\widehat{\theta }}_1\left(x\right),{\widehat{\theta }}_2\left(x\right)\right).$$

(2)

$h\left(y|\widehat{\theta }\left(x\right)\right)$ comprises both the main part and the tail part, with the latter being further divided into two segments, as illustrated in Figure 3.

4.1.1. Main part

The segment within the range of $y\le u_1\left(x\right)$ is modeled using the GMM, and its PDF is represented as $f\left(y\mid {\widehat{\theta }}_0\left(x\right)\right)$.

$$f\left(y\mid {\widehat{\theta }}_0\left(x\right)\right)=\sum _{i=1}^I{\pi }_i\left(x\right)N\left(y|{\mu }_i\left(x\right),{\sigma }_i\left(x\right)\right),$$

(3)

where

$${\widehat{\theta }}_0\left(x\right)=\left({\pi }_i\left(x\right),{\mu }_i\left(x\right),{\sigma }_i\left(x\right):i=1,2,\dots ,I\right).$$

(4)

$N\left(y|{\mu }_i\left(x\right),{\sigma }_i\left(x\right)\right)$ represents the Gaussian PDF with mean ${\mu }_i\left(x\right)$ and standard deviation ${\sigma }_i\left(x\right)$, where ${\pi }_i\left(x\right)$$\left(0<{\pi }_i\left(x\right)<1\right)$ denotes the weight. The corresponding CDF is expressed as $F\left(u_1\left(x\right)\mid {\widehat{\theta }}_0\left(x\right)\right)=P\left\{Y\le u_1\left(x\right)\mid {\widehat{\theta }}_0\left(x\right)\right\}$.

4.1.2. Tail part

For the segment within the range of $u_1\left(x\right)<y\le u_2\left(x\right)$, the PDF is a mixture of $f\left(y\mid {\widehat{\theta }}_0\left(x\right)\right)$ and $g_1\left(y\mid {\widehat{\theta }}_1\left(x\right)\right)$, as given by the second formula in (1). For the segment within the range of $y>u_2\left(x\right)$, the PDF is a mixture of $f\left(y\mid {\widehat{\theta }}_0\left(x\right)\right)$ and $g_m\left(y\mid {\widehat{\theta }}_m\left(x\right)\right)$, $m=1,2$, as described by the third formula in (1). Here, ${\alpha }_m\left(x\right)$ represents the mixture weight, and $g_m\left(y\mid {\widehat{\theta }}_m\left(x\right)\right)$ denotes the PDF of GPD, as defined in (5).

$$g_m\left(y\mid {\widehat{\theta }}_m\left(x\right)\right)=\left\{\begin{array}{cc}{\displaystyle \frac{1}{{\beta }_m\left(x\right)}}\left(1+{\displaystyle \frac{{\xi }_m\left(x\right)}{{\beta }_m\left(x\right)}}{(y-u_m\left(x\right))}^{-{\displaystyle \frac{1}{{\xi }_m\left(x\right)}}-1}\right),\hfill & {\xi }_m\left(x\right)\ne 0,\hfill \\ {\displaystyle \frac{1}{{\beta }_m\left(x\right)}}{\mathrm{e}}^{{\displaystyle \frac{y-u_m\left(x\right)}{{\beta }_m\left(x\right)}}},\hfill & {\xi }_m\left(x\right)=0,\hfill \end{array}\right.$$

(5)

where

$${\widehat{\theta }}_m\left(x\right)=\left({\beta }_m\left(x\right),{\xi }_m\left(x\right),u_m\left(x\right)\right).$$

(6)

${\beta }_m\left(x\right)>0$. $y\ge u_m\left(x\right)$ when ${\xi }_m\left(x\right)\ge 0$, and $u_m\left(x\right)\le y\le u_m\left(x\right)-{\beta }_m\left(x\right)/{\xi }_m\left(x\right)$ when ${\xi }_m\left(x\right)<0$. The corresponding CDF is expressed as $G_m\left(u_j\left(x\right)\mid {\widehat{\theta }}_m\left(x\right)\right)=P\left\{Y\le u_m\left(x\right)\mid {\widehat{\theta }}_m\left(x\right)\right\}$, $j=1,2$.

The motivations for dividing the tail part into two segments and considering the coupling relationship between different segments are as follows. The variation trend of the tail probability distribution of the measured delay indicates that the tail initially decreases slowly. However, as the delay exceeds a certain threshold, the tail probability begins to fluctuate and subsequently decreases at a faster rate. Although the GPD exhibits heavy-tailed characteristics, a single GPD is insufficient to capture the variation characteristics of the tail probability. Therefore, the tail part is divided into two segments, each of which is fitted using distinct GPDs. Additionally, according to the reference [13], fitting the probability distribution of actual data using a mixture of multiple probability distributions yields good results. Hence, the coupling relationship between the segments fitted by the GMM and GPDs is considered.

4.2. Parameter estimation for PDF

For the PDF $h\left(y\mid \widehat{\theta }\left(x\right)\right)$, the parameter vector $\widehat{\theta }\left(x\right)=\left({\alpha }_1\left(x\right),{\alpha }_2\left(x\right),{\widehat{\theta }}_0\left(x\right),{\widehat{\theta }}_1\left(x\right),{\widehat{\theta }}_2\left(x\right)\right)$ is estimated in two stages, as illustrated in Figure 4. In the first stage, the thresholds $u_1\left(x\right)$ and $u_2\left(x\right)$ in ${\widehat{\theta }}_1\left(x\right)$ and ${\widehat{\theta }}_2\left(x\right)$ are determined. In the second stage, the remaining parameters, which constitute the parameter vector ${\widehat{\theta }}^{\prime }\left(x\right)$, are estimated via MDN.

$${\widehat{\theta }}^{\prime }\left(x\right)=\left({\alpha }_1\left(x\right),{\alpha }_2\left(x\right),{\widehat{\theta }}_0\left(x\right),{\beta }_1\left(x\right),{\xi }_1\left(x\right),{\beta }_2\left(x\right),{\xi }_2\left(x\right)\right).$$

(7)

4.2.1. First stage: Determine the thresholds $u_1\left(x\right)$ and $u_2\left(x\right)$.

The threshold $u_1\left(x\right)$ is utilized to divide the main and tail parts of the delay PDF under the condition of packet length x. Based on the characteristics of heavy-tailed distributions, significant fluctuations typically occur around the inflection point between the main and tail parts. Consequently, $u_1\left(x\right)$ is defined as the delay value at which the absolute value of the differential tail probability attains its maximum. It is computed according to Algorithm 1. The inputs include the delay dataset for packets with length x, $\{(x,y_1),(x,y_2),\dots ,(x,y_{N_x})\}$, the size $\left(T+1\right)$ of the tail probability sequence, and the interval $\delta$ between adjacent delays in the tail probability sequence, such as $\delta =0.5$ ms. The outputs are the threshold $u_1\left(x\right)$ and the tail probability sequence $\{{\overline{H}}_i:i=0,1,\dots ,T\}$.

The threshold $u_2\left(x\right)$ is utilized to divide the tail part into two segments. Analysis of the measured delay dataset reveals that the tail initially decreases slowly and subsequently exhibits a faster rate of decrease. This transition occurs at a convex bend. Consequently, $u_2\left(x\right)$ is defined as the delay value corresponding to the starting point of the convex bend with the maximum vertical drop, as illustrated in Figure 5. It is determined through Algorithm 2, with inputs including T, $\delta$, $u_1\left(x\right)$, and the tail probability sequence $\{{\overline{H}}_i:i=0,1,\dots ,T\}$.

Algorithm 1: Calculate threshold $u_1\left(x\right)$

Algorithm 2 comprises three steps. In the first step, the starting points of the convex bends in the tail probability sequence $\left\{{\overline{H}}_i:i=0,1,\dots ,T\right\}$ are identified in lines 1-11 of Algorithm 2, and these points are also referred to as the starting points of convex point clusters. A convex point cluster is defined as a group of consecutive convex points. The point $\left(i,{\overline{H}}_i\right)$, $i=1,\dots ,T-1$, is termed a convex point if the following relationship (8) is satisfied.

$${\displaystyle \frac{{\overline{H}}_i-{\overline{H}}_{i+2}}{2}}+{\overline{H}}_{i+2}<{\overline{H}}_{i+1}.$$

(8)

In the second step, the convex bends in the tail probability sequence $\left\{{\overline{H}}_i:i=0,1,\dots ,T\right\}$ are identified, and their vertical drops are subsequently calculated in lines 12-26 of Algorithm 2. A set of consecutive convex points $\left\{(i,{\overline{H}}_i),(i+1,{\overline{H}}_{i+1}),\dots ,(j-1,{\overline{H}}_{j-1})\right\}$ form a convex bend if the following relationship (9) is satisfied for all $k\in \left\{i+1,i+2,\dots ,j-1\right\}$.

$${\displaystyle \frac{{\overline{H}}_i-{\overline{H}}_j}{i-j}}·\left(k-j\right)+{\overline{H}}_j<{\overline{H}}_k.$$

(9)

The vertical drop of the aforementioned convex bend is computed as ${\overline{H}}_i-{\overline{H}}_{j-1}$.

In the third step, $u_2\left(x\right)$ is determined in lines 27-33 of Algorithm 2. Specifically, $u_2\left(x\right)$ corresponds to the starting point of the convex bend with the maximum vertical drop located on the right side of $u_1\left(x\right)$.

4.2.2. Second stage: Train the remaining parameters.

The parameter vector ${\widehat{\theta }}^{\prime }\left(x\right)$ is estimated using a MDN, which comprises one input layer, four hidden layers, one output layer, and one custom layer, as illustrated in Figure 4. The custom layer is utilized to set the thresholds obtained in the first stage. The output layer consists $(6+3I)$ neurons that output the parameters in ${\widehat{\theta }}^{\prime }\left(x\right)$. The input layer contains a single neuron, into which the packet lengths $x_n$ from the dataset $\left\{(x_n,y_n):n=0,1,\dots ,N\right\}$ are fed in batches with batch size $N_b$ during each epoch. Additionally, the delays $y_n$ are used in the loss function NLPD [17], which is defined as follows.

$$\mathrm{Loss}=-{\displaystyle \frac{{\sum }_{n=k·N_b}^{(k+1)N_b-1}log\left(h\left(y_n\mid \widehat{\theta }\left(x_n\right)\right)\right)}{N_b}},$$

(10)

where $k=0,1,\dots ,$ and $(k+1)N_b-1\le N$. Moreover, the training process is carried out in multiple rounds with varying learning rates, and each round is divided into several epochs.

Algorithm 2: Calculate threshold $u_2\left(x\right)$

The aforementioned training process is adapted from reference [7]. Unlike [7], the PDF in this paper comprises a main part and a two-segment tail part. The two segments of the tail part are modeled using distinct GPDs, with consideration given to the coupling relationship between segments. The thresholds for dividing the three segments are determined based on the variation trend of tail probabilities derived from historical measured delay data, rather than through training. Meanwhile, the remaining parameters are estimated via MDN-based training.

Remark 1: The loss function NLPD is essentially a negative log-likelihood. Our objective is to minimize the negative log-likelihood, which is equivalent to performing maximum likelihood estimation (MLE). However, when the classical approach for maximizing the likelihood is applied to estimate the parameters of the delay PDF, it demonstrates poor scalability. This limitation arises from the necessity of inverting a large number of covariance matrices proportional to the number of data points [27]. In contrast, an effective alternative is the deep learning framework, which utilizes the NLPD as the loss function.

5. Experiments on WiFi Networks

In this section, we compare the proposed method TSTD with existing approaches in terms of their accuracy in estimating the delay PDF. The subsequent subsections provide a detailed description of the experimental procedures, including data preprocessing, threshold determination, model training, and performance evaluation.

5.1. Data Preprocessing

Four datasets provided by Mostafavi et al. [7] are utilized, corresponding to four types of packets with varying lengths (172 bytes, 3440 bytes, 6880 bytes, and 10320 bytes). Each dataset comprises over one million samples that record the packet length and downlink delay in software-defined radio (SDR) WiFi networks. The data collection environment for the SDR WiFi networks consists of a conference room measuring 50 m², equipped with metallic chairs, whiteboards, and screens. For collecting data related to packets with a length of 172 bytes, the coordinate of the end node location is set variably at $(1,0)$ and $(8,5)$ within a coordinate system (with scale units in meters) using the access point as the origin. This can simulate the small-scale movement of the end node. For other packet types, the end node location is consistently set at $(1,0)$. The packet generation interval was established at 10 ms. Basic information regarding the utilized datasets is presented in Table 1.

The datasets undergo preprocessing in accordance with the following steps. First, the datasets are normalized by scaling the delay values to the millisecond level. Then, standardization is performed by subtracting the mean from the normalized delay values, yielding zero-mean delay values, to avoid training errors caused by dataset biases. Finally, to enhance the convergence speed during training, the four types of packet lengths are normalized by scaling them to the range $[0,1]$.

5.2. Threshold Determination

The first stage illustrated in Figure 4 is carried out, where Algorithms 1 and 2 are employed to compute the thresholds $u_1\left(x\right)$ and $u_2\left(x\right)$ for packets with lengths $x=$ 172, 3440, 6880, and 10320 bytes, respectively. The input dataset for the algorithms takes the form $\{(x,{\dot{y}}_1),(x,{\dot{y}}_2),\dots ,(x,{\dot{y}}_{N_x})\}$, where ${\dot{y}}_n$$(n=1,2,\dots ,N_x)$ represents the delay values scaled to the millisecond level. Additionally, the other two inputs for each packet type are set as $T=200$ and $\delta =0.5$. For each packet length, the computed threshold values are shown in Table 2.

5.3. Model Training

The second stage illustrated in Figure 4 is carried out, where the parameter vector ${\widehat{\theta }}^{\prime }\left(x\right)$ is obtained by training the MDN model. In this model, the number of Gaussian distributions in (3) is set to $I=10$. The first to fourth hidden layers contain 10, 50, 50, and 40 neurons, respectively. The activation function for each neuron is set to the ’tanh’ function. For the output layer, the neurons corresponding to ${\mu }_i\left(x\right)$$(i=1,2,\dots ,10)$ have no activation function, while the neurons corresponding to ${\pi }_i\left(x\right)$$(i=1,2,\dots ,10)$ use ’softmax’ functions. The remaining neurons employ ’softplus’ functions as their activation functions. For the custom layer, the parameters $u_1\left(x\right)$ and $u_2\left(x\right)$ are fixed at the values calculated in Section 5-5.2.

The training dataset, in the form of $\{({\dot{x}}_i,{\ddot{y}}_i):i=1,2,\dots ,N\}$ is utilized, where ${\dot{x}}_i$ and ${\ddot{y}}_i$ represent the normalized packet length and delay value, respectively. This training dataset is constructed through three steps. First, the four datasets corresponding to different packet lengths are merged into a single dataset. Then, the samples in the merged dataset are randomly shuffled. Finally, N samples are randomly sampled from the shuffled dataset according to a specified sampling ratio, forming the training dataset.

The training process is carried out in four rounds with learning rates of ${10}^{-2}$, ${10}^{-3}$, ${10}^{-4}$, and ${10}^{-5}$, respectively. Each round consists of 200 epochs. In each epoch, the batch size $N_b$ is set to $1/8$ of the number of samples in the training dataset.

5.4. Performance Evaluation

The performance of TSTD is evaluated in terms of the means and fluctuation amplitudes of the delay tail probabilities, and it is compared with two existing methods. One method is GMM [15], which is one of the state-of-the-art approaches for probability density estimation. The other method is GMEVM [7], which integrates GPD with GMM.

First, Q training datasets are independently sampled with a sampling ratio of $0.8\%$, and these datasets are used to train and obtain Q PDFs. In the experiment, Q was set to 9. The mean and fluctuation amplitude of Q tail probability distributions are plotted in Figure 6, where the delay value serves as the abscissa and the common logarithm of the tail probability serves as the ordinate. The solid line represents the mean, while the shaded area represents the fluctuation amplitude. It can be observed that both GMM and GMEVM exhibit larger deviations from the actual measurements (denoted as MEAS), whereas TSTD fluctuates around the actual measurements.

Then, similar to the above case, the experiment results for the sampling ratios $6.25\%$ and $21\%$ are presented in Figure 7 and Figure 8, respectively. It can be observed that the fluctuation amplitudes of both GMEVM and TSTD decrease as the sampling ratio of the training dataset increases. TSTD always follows the actual measurements as the delay increases, while GMEVM initially exhibits heavy-tail characteristics and gradually deviates from the actual measurements. Because a single GPD is insufficient for accurately capturing the probabilistic characteristics when the tail exhibits heavy-tailed behavior accompanied by fluctuations. Even though the sampling ratio increases to $21\%$, GMM continues to exhibit significant fluctuations relative to the actual measurements. This is due to the fact that as the sampling ratio increases, the influence of tail events becomes more pronounced. Fitting heavy-tailed distributions using GMM introduces substantial errors given the inherent limitations of GMM. This issue persists even when the number of GMM components is increased, as shown in Figure 9, where the number of components is set to 20.

Next, the model is trained using the complete dataset. As shown in Figure 10, the improvement in estimation performance is not significant. In dynamic wireless environments, the model needs to be updated periodically. Introducing a data sampling strategy can effectively reduce computational costs. Thus, achieving an appropriate balance between the sampling ratio (or training time) and model performance is essential for efficient model updates. Table 3 provides a detailed summary of the average training time for Q PDFs under different sampling ratios. Adding one segment means that the model needs to train additional parameters, which also increases the training time to some extent. TSTD requires longer training time compared to GMM and GMEVM. However, the difference is not significant compared to GMEVM, especially when the size of the training dataset is larger.

Furthermore, the mean absolute error ($\widehat{MAE}$) and average fluctuation amplitude ($\widehat{AFA}$) of tail probabilities across all considered delays are utilized to evaluate the performance of estimation methods.

$$\widehat{MAE}\left(x\right)={\displaystyle \frac{\sum _{j=0}^{T}MAE\left(y=j·\delta ,x\right)}{T+1}},$$

(11)

where $MAE\left(y=j·\delta ,x\right)$ represents the mean absolute error of tail probabilities at delay y, conditioned on packet length x, considering the Q PDFs estimated from Q training datasets.

$$MAE\left(y=j·\delta ,x\right)=|{\displaystyle \frac{\sum _{i=1}^Q{\overline{H}}^i\left(y|\widehat{\theta }\left(x\right)\right)}{Q}}-MEAS\left(y|x\right)|.$$

(12)

$Q=9$. ${\overline{H}}^i\left(y|\widehat{\theta }\left(x\right)\right)$ represents the tail probability at delay y given packet length x, estimated independently using the i-th training dataset, where $i=1,2,\dots ,Q$. Additionally, $MEAS\left(y|x\right)$ denotes the measured tail probability at delay y under the condition of packet length x.

$$\widehat{AFA}\left(x\right)={\displaystyle \frac{\sum _{j=0}^{T}AFA\left(y=j·\delta ,x\right)}{T+1}},$$

(13)

where $AFA\left(y=j·\delta ,x\right)$ represents the average fluctuation amplitude of tail probabilities at delay y given packet length x, considering Q PDFs estimated from Q independent training datasets.

$$\begin{array}{cc}\hfill & AFA\left(y=j·\delta ,x\right)\hfill \\ \hfill & =|MAX\left\{{\overline{H}}^i\left(y|\widehat{\theta }\left(x\right)\right):i=1,2,\dots ,Q\right\}\hfill \\ \hfill & -MIN\left\{{\overline{H}}^i\left(y|\widehat{\theta }\left(x\right)\right):i=1,2,\dots ,Q\right\}|.\hfill \end{array}$$

(14)

$MAX\{·\}$ and $MIN\{·\}$ denote the functions used to select the maximum and minimum values, respectively, from a given set.

As shown in Table 4, for TSTD with a sampling ratio of $21\%$, both $\widehat{MAE}$ and $\widehat{AFA}$ decrease compared to GMM and GMEVM.

For some delay values, the MAE and AFA of TSTD either remain at the same order of magnitude or decrease by at least one order of magnitude compared with GMM and GMEVM. For example, the MAEs and AFAs for certain delay values are listed in Table 5 and Table 6. However, for specific delay values, the MAE of TSTD increases by one order of magnitude compared to GMM, while the AFA decreases by one order of magnitude, with the fluctuation interval still lying within that of GMM. The PDFs constructed in GMM and TSTD are used to approximate the PDF of E2E delays. In this approximation process, the loss function is defined as the mean of negative log probability densities over multiple delays, which measures the average closeness between the estimated result and the target. Therefore, compared with GMM, it is possible for TSTD to exhibit poor fitting results on a small number of points, but its overall performance remains superior.

Finally, the potential inaccuracies or noise in delay measurements are considered as factors that could influence the estimation results. To evaluate the performance of estimation methods under such conditions, random Gaussian noises with variances of 1 ms and 3 ms are respectively introduced into the delay samples. As shown in Figure 11 and Figure 12, in the low-delay domain, as the noise variance increases, the mean of each method gradually deviates from the actual measurements, while their fluctuation amplitudes remain largely unaffected. In the high-delay domain, as the noise variance increases, the mean and fluctuation amplitude of TSTD remain almost unchanged. For GMEVM, its fluctuation amplitude initially increases and then decreases, while its mean remains nearly constant. For GMM, both its mean and fluctuation amplitude are relatively stable, with the latter continuing to decrease slightly. This is because a higher proportion of data lies in the low-delay domain, where changes in these data can significantly affect the mean of the fitting results. Additionally, the added noise tends to smooth out sharp fluctuations, ultimately reducing the fluctuation amplitudes. In summary, compared with GMM and GMEVM, the multiple-segmentation feature of TSTD makes it more robust when measurement data are slightly inaccurate or contain low levels of noise.

6. Conclusion and Future Work

In this paper, a PDF estimation method for the E2E delay in WiFi networks is proposed. The delay PDF is divided into three segments: the first segment, fitted by GMM, constitutes the main part of the PDF, while the last two segments, fitted by two different GPDs, constitute the tail part of the PDF. The thresholds for dividing the three segments are calculated first, and the other parameters of the PDF are subsequently obtained through training using MDN. Experimental validation demonstrates that this approach achieves results closer to measurement data compared to existing methods. Specifically, the mean absolute errors and average fluctuation amplitudes of tail probabilities at certain delay values decrease by at least one order of magnitude. Moreover, the multiple-segmentation feature of the proposed method enhances its robustness in situations where measurement data are affected by low levels of Gaussian noise.

In the future, further refinements will focus on integrating additional network parameters as conditioning variables to enhance both the model’s accuracy and adaptability. Potential parameters encompass modulation and coding scheme (MCS), packet arrival intervals, signal-to-noise ratio (SNR), and packet loss rate. Moreover, assessing the model’s generalizability beyond WiFi networks, such as in 5G and IoT networks, represents a critical research direction. Ultimately, the estimation results will provide a robust foundation for devising deterministic delay control strategies.