Introduction

During researching and modeling information processes in telecommunication systems and computer networks, it is necessary to take into account the fractal properties of telecommunication traffic. In particular, taking into account fractal properties significantly affects the results of simulation modeling of information exchange processes, and, accordingly, it allows predict the development of events in information processes that are associated not only with the movement of information in telecommunications systems [1][2][3][4][5][6]. All these sources are related to the problem of determining the fractal dimension of time series, and the dates of publications vary from 1994 to 2021, which indicates that the problem of rapid and accurate determination of fractal dimensions is not finally solved and needs further research.

Often, to determine the fractal characteristics of time series, used methods are based on the simulation of random walks [7][8][9][10][11], which depends on the standard deviation of the random variable and the length of the series. In particular, such an indicator is R/S analysis [12][13][14][15][16]. It is believed that R/S analysis for searching the self-similarity has a higher accuracy for long series with large cumulative sums [6]. Also, it is noted that determining the self-similarity coefficient depends on the length of the researched series [17]. It is not always possible to obtain long time series for its analysis. For example, time series properties change before it will be able to be used the calculated coefficient, or the process itself is too short to obtain enough data for reliable analysis. Therefore, it can be seen that the flow of scientific works on methods for obtaining the self-similarity coefficient of the time series does not decrease. Eventually it demonstrates the relevance of the given problem.

Related Works and Problem Statement

So, during solving problems of forecasting and optimizing the operation of telecommunication networks, it is almost necessary to take into account the fractal nature of traffic on the Internet [18][19][20][21].

The known criteria by which fractal dimension is determined, for example, the Hurst exponent [22], have significant errors and deviations for individual implementations, so it is advisable to develop new methods for assessing the fractal characteristics of the signal under study, such researches are observed in a number of publications, for example, in [23] was proposed method using wavelet analysis. Some review publications, which consider and compare several methods, use the definition of the fractal dimension of individual short-term implementations, based on which determine the properties of this signal [24]. But there are cases of ergodic signal with a known theoretical justification of its probability density distribution function. In such cases, short implementations have a high probability of deviations, but such implementations allow us to estimate the weights of the known probability distribution. Exactly in such cases it is useful to have high-quality means of determining the fractal dimension based on the probability distribution function, which is the reason for the constant search and development of new better methods for determining the fractal dimension. In the processes of information transmission in telecommunication systems and networks, the search for the fractal dimension of network traffic is widely used. To do this, use different methods and approaches aimed at reducing random deviations when performing calculations on relatively small implementations of numerical series [19-21, 23, 24]. Given the lack of input data, additional data on the nature of the studied time series can significantly improve the situation. For example, with the help of theoretical justifications for modeling time series, they can be realized on the basis of random processes with given probability distributions, such as Poisson process [25, p. 32], Markov chains [25, p. 89], queues based on Pareto process [26] and others. Each of the distributions has its own scope and is based on accepted hypotheses and experimental testing on long-term implementations. However, real processes have only approximate distributions to the theoretical ones, it is worth remembering that the theoretical flow of data per unit time is limited by the bandwidth of the input channel of the network. Therefore, the approximation can be improved by taking into account the real probability distributions obtained experimentally.

The authors proposed a method for determining the fractal dimension of a numerical sequence by the probability distribution of the values of its elements (Probabilistic method) in [27,28]. In the mentioned work, the problem of obtaining a mathematical relation to obtain the expected value of the fractal dimension of network traffic based on the knowledge of the probability density distribution p(x) was set and solved. The results can be used to predict changes in network traffic, as well as to improve the modeling of fractal random sequences, queuing systems [29], simulating the traffic of telecommunications networks [30].

The purpose of this work is an experimental study of the quality of the proposed method for determining the fractal dimension of network traffic by its probabilistic properties.

A method for determining the fractal dimension of network traffic by its probabilised properties 2.1. Theoretical substantiation and essence of the method

Modern traffic management systems take into account its fractal properties. Network traffic can be represented as a time series. In systems that use packet information exchange, it is more convenient to use a binary time series, and the representation of traffic at the level of "data packet is present", "there is no data packet". In this work, a probabilistic analysis of such a time series is performed in order to highlight the probabilistic characteristics that reveal not only the intensity of such traffic, but also its fractal properties.

The following relations are valid for the probability density distribution p(x):

∫ 𝑝(𝑥)𝑑𝑥 +∞ −∞ = 1, 𝑝(𝑥) ≥ 0, 𝑀(𝑥) = ∫ 𝑥 • 𝑝(𝑥)𝑑𝑥 +∞ −∞(1)

where M(x)is the mathematical expectation of the random variable x. A numerical sequence consists of a series of implementations xi, where iis the serial number of the element. By definition, the Minkowski dimension is the value of the following limit:

lim 𝜀→0 ln(𝑁 𝜀 ) ln(𝜀) ,

where εis the diameter of the coating element; Nεis their number. In practice, the problem is solved geometrically by covering the studied figure with squares (cubes), where its side is taken as the diameter. In [1,31] there is a variant in which the coverage is replaced by rectangles of width ε and height, which is minimal to cover the area of the graphical representation of the numerical series. When using rectangles, the number of figures of the coating is replaced by the coverage area:

𝑆(𝜀) = ∑ ℎ 𝑘 𝑁 𝜀 𝑘=0 ,(2)

where hkis the height of the corresponding rectangle k.

If take ε > 1 for the number of discrete samples of the numerical series xi, then hk is searched according to the following algorithm:

ℎ 𝑘 = max(𝑥 𝑖 ) − min (𝑥 𝑖 ) ,

where 𝑖 = 𝑘𝜀, … , (𝑘 + 1)𝜀 − 1.

For an inclined straight line, the coating with a decrease ε in 2 times reduces the area also twice. For a figure that covers a plane (for example, the Hilbert or Piano curve), when ε changes, the value of the area does not change and covers the entire part of the plane. This corresponds to the fractal dimension of the plane D = 2. Therefore, the fractal dimension is expressed through the area of coverage [2,3] by relation (4), if the width of the rectangle is changed γ times:

𝐷 = 2 − log 𝛾 (𝑆(𝜀 • 𝛾)/𝑆(𝜀)).(4)

For a significant number of samples of the numerical sequence N, it is possible to replace the values of the heights of the rectangles of coverage hk by their mathematical expectation M(hk). But the mathematical expectation of the height of the rectangle depends on the number of samples per partition ε, and this dependence is nonlinear. Therefore, we denote the mathematical expectation of the height of a rectangle depending on its width ε as М(ε), where the mathematical expectation of the area of an individual rectangle will be expressed as the product of its width and the expected height (5):

𝑆 𝑘 (𝜀) = 𝑀(𝜀) • 𝜀.

(

)5

Then have the expression for determining the fractal dimension through the coverage of rectangles of formula (2), taking into account (5):

𝑆(𝜀) = ∑ 𝑆 𝑘 (𝜀) 𝑁 𝜀 𝑘=0 ⇒ 𝑆(𝜀) = ∑ 𝑀(𝜀) • 𝜀 𝑁 𝜀 𝑘=0 ⇒ 𝑆(𝜀) = 𝑀(𝜀) • 𝜀𝑁 𝜀 ⇒ 𝑆(𝜀) = 𝑀(𝜀) • 𝑁, (6)

where N is the total number of elements of the series, which is divisible by ε (although for N >> ε this requirement can be neglected and incomplete last rectangles of coverage can be accepted).

If use the obtained relation (6) to the dimension value, obtain:

𝐷 = 2 − log 𝛾 ( 𝑀(𝜀 • 𝛾) 𝑀(𝜀) ),(7)

in which the number of numbers in the implementation of series N has been reduced and is no longer used. Unfortunately, the mathematical expectation of the height of the rectangle (3) is not known. Therefore, we obtain from the mathematical expectation of the quantity itself. To do this, find the probability q3(h) of the height of the rectangle for ε = 2:

𝑞 2 (ℎ) = ∫ 2𝑝(𝑥) • 𝑝(𝑥 + ℎ)𝑑𝑥 +∞ −∞ ,(8)

an additional factor of 2 is used here, because the resulting probability turns out in two independent cases when the smaller number is obtained earlier and later.

The general expression for q3(h) is obtained if all subsequent realizations of the number x, by the number ε-2, will be located on the interval [x, x + h]. To do this, the probability of ∫ 𝑝(𝑥)𝑑𝑥 𝑥+ℎ 𝑥 falling into this interval must be increased to the power of ε-2 and multiplied by their possible combinations of the location of the limit values x and x + h, as which acts binomial coefficient 𝐶 𝜀 𝜀−2 = 𝜀(𝜀 − 1)/2:

𝑝 𝜀−2 (𝑥, 𝑥 + ℎ) = 𝜀(𝜀 − 1) ( ∫ 𝑝(𝑡)𝑑𝑡 𝑥+ℎ 𝑥 ) 𝜀−2 ,(9)

Therefore, for q3(h), the expression is written with (8) taking into account (9) as the following integral (10):

𝑞 3 (ℎ) = 𝜀(𝜀 − 1) ∫ 𝑝(𝑥) • 𝑝(𝑥 + ℎ) ( ∫ 𝑝(𝑡)𝑑𝑡 𝑥+ℎ 𝑥 ) 𝑑𝑥 +∞ −∞ .(10)

Formula ( 10) in fact provides information about the probability density distribution of the height of the rectangle of coverage with a width of ε samples. The value of the expression is not valid for ε < 2, because for a single width of the rectangle, it should cover only one point, so its height will always be zero. In other cases, zero or less samples cannot be taken, because the number of samples ε is a natural number.

Now that there is a value for the probability distribution of the width of the rectangle (10), we can write the value of the mathematical expectation of this height (11):

𝑀(𝜀) = 𝜀(𝜀 − 1) • ∫ ℎ ( ∫ 𝑝(𝑥) • 𝑝(𝑥 + ℎ) ( ∫ 𝑝(𝑡)𝑑𝑡 𝑥+ℎ 𝑥 ) 𝜀−2 𝑑𝑥 +∞ −∞ ) +∞ 0 𝑑ℎ. (11)

For the mathematical expectation M(ε) the integration limit is taken on the interval [0; +∞], because the height of the rectangle cannot have a negative value.

Finally, using the obtained mathematical expectation of the height of the rectangle (11) for expression (7), we have an estimate of the fractal dimension based on the density of the distribution (1):

𝐷 = 2 − log 𝛾 ( 𝛾(𝜀 • 𝛾 − 1) ∫ ℎ (∫ 𝑝(𝑥) • 𝑝(𝑥 + ℎ) (∫ 𝑝(𝑡)𝑑𝑡 𝑥+ℎ 𝑥 ) 𝜀•𝛾−2 𝑑𝑥 +∞ −∞ ) +∞ 0 𝑑ℎ (𝜀 − 1) • ∫ ℎ (∫ 𝑝(𝑥) • 𝑝(𝑥 + ℎ) (∫ 𝑝(𝑡)𝑑𝑡 𝑥+ℎ 𝑥 ) 𝜀−2 𝑑𝑥 +∞ −∞ ) +∞ 0 𝑑ℎ ). (12)

Sometimes it is advantageous to use a non-multiple partition relation, so the following notation of expression (12) will be more useful:

𝐷 = 2 − log 𝜀 2 /𝜀 1 ( 𝜀 2 (𝜀 2 • 𝛾 − 1) ∫ ℎ (∫ 𝑝(𝑥) • 𝑝(𝑥 + ℎ) (∫ 𝑝(𝑡)𝑑𝑡 𝑥+ℎ 𝑥 ) 𝜀 2 −2 𝑑𝑥 +∞ −∞ ) +∞ 0 𝑑ℎ 𝜀 1 (𝜀 1 − 1) • ∫ ℎ (∫ 𝑝(𝑥) • 𝑝(𝑥 + ℎ) (∫ 𝑝(𝑡)𝑑𝑡 𝑥+ℎ 𝑥 ) 𝜀 1 −2 𝑑𝑥 +∞ −∞ ) +∞ 0 𝑑ℎ ).(13)

A formula for obtaining a mathematical expectation of the fractal dimension of a sequence of random numbers with a known probability density distribution ( 13) is obtained.

The obtained formula ( 13) is transient for estimating the fractal dimension of a binary series, which consists only of elements that are equal to 0 or 1 (Fig. 1).

Figure 1: Visualization of zero coverings on a discrete binary number series

To estimate the fractal dimension of such a binary series, need to determine the probability of obtaining a sequence of a given length ε zeros p0(ε) and units p1(ε):

𝑝 0 (𝜀) = λ 0 𝜀 , 𝑝 1 (𝜀) = λ 3 𝜀 ,

where λ0 is the probability of maintaining the zero state, λ3 is the probability of maintaining the single state. According to the defined probabilities, the mathematical expectation of the degree of coverage n elements of the sequence will be expressed by the following formula:

𝑀(𝑛) = 1 − (1 − λ 3 )λ 0 𝑛 + (1 − λ 0 )λ 3 𝑛 (1 − λ 0 ) + (1 − λ 3 ) .

If, for convenience, go to probability designations to change the current state in the series from zero to one and vice versa, λ1=1-λ0 and λ2=1-λ3, then in (13) can use discrete sums instead of integrals. Determining the mathematical expectation of the degree of coverage M(n), allows to use the following formula to determine the fractal dimension of the binary sequence:

𝐷(𝑛 1 , 𝑛 2 , λ 1 , λ 2 ) = 1 + 𝑙𝑛 ( λ 1 + λ 2 − λ 2 (1 − λ 1 ) 𝑛 1 − λ 1 (1 − λ 2 ) 𝑛 1 λ 1 + λ 2 − λ 2 (1 − λ 1 ) 𝑛 2 − λ 1 (1 − λ 2 ) 𝑛 2 ) 𝑙𝑛(𝑛 1 /𝑛 2 ) ⁄ .

The fractal dimension search expression contains not only the transition probabilities but also the coverage lengths. This leads to ambiguity in the resulting dimension. To get rid of the ambiguity, direct n1 to n2 and then direct the number of elements to 1. In result of finding the boundaries, we obtain the final formula for determining the fractal dimension of a binary sequence by transition probabilities [27]:

𝐷(λ 1 , λ 2 ) = 1 − λ 2 (1 − λ 1 )𝑙𝑛(1 − λ 1 ) + λ 1 (1 − λ 2 )𝑙𝑛(1 − λ 2 ) 2λ 1 λ 2 . (14)

Due to the fact that the binary time series is given by the probabilities of changing or leaving the state, it was possible to express the fractal dimension analytically (14). Visually, this can be assessed in Fig. 2. Unfortunately, the coating method has an effect on the value of the fractal dimension, so to compare the results you need to make the transition from the resulting dimension to others by regression.

As a result, we obtained a method for determining the fractal dimension of network traffic by its probabilistic properties: Stage 1. The binary time series for which it is necessary to define fractal properties accumulates. Stage 2. The probabilistic properties of the series are evaluated. The accumulated data calculates the number of occurrences of pairs of consecutive elements n00 : "0, 0" and n01 : "0, 1". We obtain λ1 as the probability of change of state from "0" to "1" by the definition of probability λ1 = n01 / (n01 + n00).

Similarly, the probability of changing the state from "1" to "0" is estimated as λ2 = n10 / (n10 + n11), where n10 is the number of pairs of consecutive elements "1, 0"; n11is the number of pairs of consecutive elements "1, 1".

Stage 3. The fractal dimension of a binary sequence is determined by formula (14).

Figure 2: Dependence of the fractal dimension of a binary series on the probabilities of changes in the state of the system For practical use of the results there is a need to create a binary sequence generator with specified probability characteristics. This question is discussed below.

Experimental study of the quality of the method

To conduct experiments, it was decided to simulate network traffic with predefined properties. It was decided to model based on the theory of Markov processes, which is often used to model the traffic of different queuing systems [32][33][34][35][36][37][38]. For generation fractal binary traffic used Markov chain, shown in Fig. 3.

In this work, a binary time series was created to simulate network traffic, the persistence of which is regulated by setting the probabilities of state change to the opposite λ1, λ2 (Fig. 3). 0 1

λ 2 λ 1 λ 0 λ 3 p 0 p 1

Figure 3: Markov chain for generating fractal binary traffic

This generator is characterized by states 0 or 1, and the probabilities of being in these states as p0 = λ2/(λ1 + λ2) and p1 = λ2/(λ1 + λ2), where λiis the probabilities of the corresponding transitions [27]. The traffic intensity of such a generator will be within [0, 1] and will be equal to the probability of obtaining the output of the generator 1: p1. The algorithm of operation of such a generator is shown in Fig. 4. Generating traffic with intensity ½ begins with setting the probability of maintaining the state λ. The algorithm contains a variable to save the previous state. The cycle repeats the generation of a pseudorandom number from the range [0; 1) with a uniform distribution, for which a comparison is made with a given probability λ. When passing the test for comparison, the state is preserved, and the output is given the value of the previous state, otherwise the state changes to the opposite. It is expected that the generated binary traffic according to this algorithm has a controlled fractal dimension according to relation (14).

In order to determine the fractal properties of the binary series obtained by the generator from Figs. 3, an experimental measurement of the Hurst exponent by R/S analysis was performed, the results of which are shown in Fig. 5.

Fig. 5 contains the results of the analysis of binary traffic with intensity p1 = 0.5. To do this, it suffices to fulfill the condition of equality of transition probabilities λ = λ0 = λ3, or equivalent to λ1 = λ2. This indicator is responsible for the persistence of the time series and affects its fractal properties. Also, to calculate the Hurst exponent H, need to select a series of random wanderingscumulative sums. The number of steps of the cumulative sums is shown by a separate axis L. The graph shows that as the length of the cumulative sums decreases, the graph goes to a straight line that connects the unit and zero values of the Hurst exponent. Сonversely, if the cumulative sums are long enough, the Hurst exponent goes to a value of 0.5 without reflecting the persistence of the time series.

The Hurst exponent was calculated according to the R/S analysis algorithm, which is implemented by the following program code [31]: The algorithm used for the experiment is shown in Fig. 6. The algorithm takes at the input the sequence to be analyzed, and lagsis the value of shifts in the calculation of cumulative sums. For each point of the graph, which corresponds to the probability of maintaining the state and range of lengths of cumulative sums, 100 times a random series of 10000 values is generated, for which Hurst exponents are calculated with the specified lengths of cumulative sums. The graph shows the average value for these 100 measurements. Boundary results need better clarity. In Fig. 7 shows the result of using the minimum cumulative sums in 1 and 3 consecutive values: As a result, we have a curve that is an approximation of the linear relationship between the probability of maintaining the previous state λ and the Hurst exponent. Also, the markup of deviations shows a small discrepancy between the results from experiment to experiment. However, this is the result of a deliberately election long series of experiments in order to obtain reliable experimental results.

1 # Input data: 2 #

The following Fig. 8 is a representation of another extreme position (in Fig. 5, this trend is developing "deep"), when the cumulative sums are long: Here it is obvious that there is a significant deviation in some measurements, which affected the bandwidth of 2σ. This indicates that the measurement of the Hurst exponent is extremely inaccurate. In this case, in a wide range of probabilities of preservation of the previous state, about [0.1; 0.9], R/S analysis showed that the series is random.

Comparison of graphs in Fig. 7 and Fig. 8, as a result of the analysis of the results of Fig. 5, leads to the conclusion that the result depends on the user's choice of parameters, in particular the length of the series or the lengths of the cumulative sums.

It is proposed to use the Probabilistic approach, which is based on direct measurement of the probability of state change from 0 to 1 (r0) and from state 1 to state 0 (r1) according to the formula for determining fractal dimension (14), which is implemented in Python 3 programming language: where, series is a one-dimensional array containing a binary time series. Using this code to determine the fractal properties of the time series is shown in Fig. 9.

1 r0 = (series[0:-2] < series[1:-1]).

The results of the research of the quality of the Probabilistic method and its comparison with the well-known method of R/S analysis are shown in Table 1.

To construct each row of the table, 100 generations of a binary sequence with a length of 10000 elements with the corresponding value of λ were performed. For each of them, the Hurst exponent H was determined by R/S analysis and by based on an estimate of the probability of transition from 1 to 0 and from 0 to 1 (Probabilistic method, PM).

Thanks to hundreds of independent tests, it is possible to determine the average Hurst exponent (15) standard deviation σ (16) and percentage of deviation (17): 0,06 0,115130822 0,030678064 0,005365203 0,001102508 4,660093493 3,593801913 3 0,11 0,197424495 0,057299125 0,006720639 0,001636951 3,404157102 2,856851724 4 0,16 0,266660733 0,084636434 0,008166731 0,001946542 3,062592345 2,299887334 5 0,21 0,321224610 0,113072840 0,009107072 0,002310433 2,835110483 2,043314238 6 0,26 0,367964061 0,143068708 0,009842329 0,002886502 2,674807286 2,017564092 7 0,31 0,404989816 0,173980449 0,010944700 0,003172788 2,702463230 1,823646757 8 0,36 0,435730431 0,206981323 0,011123077 0,003080814 2,552742787 1,488450246 9 0,41 0,462658106 0,240125774 0,011055708 0,003745437 2,389606470 1,559781461 10 0,46 0,485639655 0,277501570 0,009411981 0,003754746 1,938058694 1,353054180 11 0,51 0,502866839 0,313936619 0,011826159 0,004215003 2,351747645 1,342628967 12 0,56 0,520859605 0,354596924 0,011126878 0,004295506 2,136252974 1,211377335 13 0,61 0,540472278 0,398324387 0,010182414 0,004817891 1,883984610 1,209539751 14 0,66 0,561685180 0,444444295 0,012015554 0,004776983 2,139197425 1,074821662 15 0,71 0,587349580 0,494695839 0,010646189 0,004218451 1,812581414 0,852736347 16 0,76 0,619086392 0,549212501 0,010121335 0,005005190 1,634882623 0,911339444 17 0,81 0,667516915 0,609945932 0,010737486 0,005076275 1,608571380 0,832250120 18 0,86 0,737425650 0,680389339 0,009941918 0,004865064 1,348192641 0,715041204 19 0,91 0,836641573 0,762035440 0,010570512 0,004452206 1,263445799 0,584251869 20 0,96 0,949146972 0,866270314 0,007584074 0,005016391 0,799041091 0,579079216 Averages: 0,009416 0,003541 2,597437 1,852847

To build each row of the table, 100 binary sequence generations with a length of 10000 elements with an appropriate value of λ were performed. There are σis the standard deviation σ (16) and Deviationsis the percentage of deviation (17).

The constructed table allows to estimate the accuracy of the found Hurst exponents by last two columns (Table 1). According to these columns, the relative deviation for the Probabilistic method is smaller. This is because R/S analysis requires to estimate the "scatter" of R, which is the cumulative sum of several elements of the studied sequence. It is clear that such a value requires more measurements of traffic values and requires estimating the values with the appropriate accuracy, i.e. more experiments.

𝐻= 1 100 ∑ 𝐻 𝑖 , 100 𝑖=1(15)𝜎 = √ 1 100 ∑(𝐻 𝑖 − 𝐻) 2 , 100 𝑖=1(16)Deviations = 𝜎 𝐻 • 100%.(17)

R/S analysis gives different values for different lags of the cumulative sum and the larger the step of the cumulative sum, the more different the value of the Hurst exponents determined by it. This indicates the ambiguity of the method and the possibility of manipulating the results by choosing the size of the cumulative sums. The proposed Probabilistic method does not have this disadvantage. The developed method is unambiguous and does not depend on the period of consideration, because it does not use cumulative sums.

Also, as can be seen from Table 1, developed Probabilistic method has a lower percentage deviation from the mean value of the Hurst exponents, namely 1.8% as opposed to the R/S analysis, which has a deviation 2.5%. Therefore, the proposed method has a greater accuracy in determining the fractal dimension.

Conclusions

This experimental study investigates the quality of work of the previously proposed method of determining the fractal dimension of network traffic by its probabilistic properties. Also, the results of the experimental comparison of the work quality of the proposed Probabilistic method compared to the method based on R/S analysis.

A simulated binary time series using Markov process theory was used to conduct the experiments. To simulate network traffic a binary time series was created with persistence regulated by setting up the probabilities of changing states "0" and "1" to the opposite ones.

As shown by experiments on simulated network traffic, the proposed Probabilistic method has greater accuracy and unambiguity of results regardless of the length of the studied series in contrast to R/S analysis. R/S analysis gives different values for different sizes of the cumulative and the larger the step of the cumulative sum, the more diverges the value of the Hurst exponent determined by it. This indicates the ambiguity of the method and the possibility of manipulating the results by choosing the different lags of the cumulative sums. The proposed Probabilistic method does not have this disadvantage. The developed method is unambiguous and does not depend on the period of consideration, because it does not use cumulative sums.

Also, the developed method has a lower percentage of deviation from the mean value of the Hurst exponent and therefore has higher accuracy in determining the fractal dimension -R/S analysis has a deviation of 2.5%, and the developed method 1.8%.