This paper proposes an approach for the estimation of a time-varying Hurst exponent to allow accurate identification of multifractional Brownian motion (MFBM). The contribution provides a prescription for how to deal with the MFBM measurement data to solve regression and classification problems. Theoretical studies are supplemented with computer simulations and real-world examples. Those prove that the procedure proposed in this paper outperforms the best-in-class algorithm.

Time-changing fractional dynamics are identified in signals measured in systems of various nature, e.g., biological, technical, financial, etc. One of the models used to describe this phenomenon is the multifractional Brownian motion (MFBM), considered a generalization of the classical fractional Brownian motion (FBM). The MFBM is characterized by a time-dependent Hurst exponent [$H(t)$], in contrast to the FBM where $H$ is constant. In this article, an artificial intelligence-based method dedicated to the estimation of time-varying $H(t)$ from experimental data has been introduced. The algorithm uses a deep convolutional neural network (CNN) with an encoder–decoder architecture, the autocovariance function (ACVF), and the Euclidean distance matrix of the considered process. We demonstrate how the procedure can be applied when noisy and only a limited amount of measurement information is available for estimation. Simulation studies with different time-changing $H(t)$ show that the proposed algorithm enables estimation of the Hurst coefficient with improved accuracy in relation to Coeurjolly’s method, currently perceived as the best in class. We propose a discriminating procedure that can be used for the classification of trajectories for cases with constant and time-varying fractional dynamics, especially important for context labeling of experimental data. The simulation study is supplemented by experimental studies of signals recorded in a hydrogel system. Although the proposed methodology is demonstrated for the MFBM, the designed and validated approach is of universal meaning and can be applied for systems of any nature where the time-changing fractional dynamics govern a temporal evolution.

## I. INTRODUCTION

Complexity identified in measured signals reflects the complex interrelations (structural and functional) valid for observed systems.^{1,2} This property does not depend on the scale and nature of the system; i.e., one can observe complex behavior in geological, meteorological, financial, biological, and engineering objects perceived in a macro-world context, e.g., Refs. ^{3–8}. Nevertheless, the complex patterns are also visible when we characterize the micro-to-nano world, e.g., plasma processes, genomics, and materials science.^{9–11} Fractional calculus is one of the mathematical approaches used to model this class of systems and processes, but applications in measurement science and engineering benefit from simple and efficient procedures for the depiction of the rules and properties governing complexity. Fractional Brownian motion (FBM) with Hurst exponent $H$ is a stochastic process that can entail long-ranged temporal correlations^{12–15} and has been broadly employed across multiple disciplines.^{16,17} It has been established experimentally and defined theoretically to depict complex interrelations and system dynamics.^{18,19} Typically, the estimation of the constant $H$ on a given scale or within a multiscale view is performed to reveal and label the system state, e.g., healthy vs asthmatic.^{4,20} However, real systems often evolve in time such that the Hurst index is not constant along the whole time series. For example, persistent or anti-persistent processes, where the state changes with time, have been observed in systems as diverse as the motion of proteins in living cells,^{21,22} nanoscale particles in the cytoplasm,^{23,24} and the degradation of Li-ion batteries.^{25} It has been shown that the index $H$ can change with time, and thus, a more advanced description [with a reconstructed variation in $H(t)$] should be applied. An example of such a process is the generalization of FBM and the multifractional Brownian motion (MFBM),^{26–29} as a theoretical model, which is valid for the description of real-world systems of any nature and across various scales. A mathematical consequence of relaxing the stationarity constraint in the MFBM model lies in the ability to control the local irregularity at the output, which is valid for physical systems. Thus, the MFBM process provides a useful description for a host of continuous and non-stationary natural signals, exhibiting physical complexity encoded in observed systems and, in turn, broadening the field of its application. One can imagine a pointwise irregularity evoked by erosion phenomena in pricing, mountain rock, glacier degradation, etc.^{30–34} Some of these components cannot be directly measured, and they need to be reconstructed from experimental data; e.g., a projection of the relation between cell’s configuration and performance in a wireless communication network needs to assume, among other things, time-varying environmental changes, which profile the output of a technical system operating in a given geographical location.^{35} Furthermore, a challenge in signal processing is to reliably estimate a varying $H(t)$ as the fingerprint of the dynamics when the system is corrupted by a noisy environment using a limited amount of *a priori* information. In this article, a prescription for solving this challenge is provided for the MFBM, considered the classical model with time-changing fractional dynamics. In contrast to FBM, the MFBM assumes a time-dependent Hurst exponent. In the literature, one can find various algorithms for the estimation of a time-dependent Hurst exponent.^{31,36–41}

We develop an approach based on artificial intelligence (AI) that is verified with numerical simulations and a real-world system: signals recorded in a microscale experiment with polystyrene nanobeads placed in an agarose hydrogel. The proposed technique utilizes deep convolutional neural networks (CNNs) with encoder–decoder architecture.^{42–44} The presented results, i.e., the method of time-varying $H(t)$ estimation, and qualitative and quantitative observations performed for a theoretical and real-world case can be translated to the real-world system/process of any nature and scale where the fractional dynamics can be explained by the time-changing Hurst exponent.

The application of the AI and machine learning (ML) techniques to signal processing brings vivid interest of researchers and engineers^{45,46} since this approach is efficient to the exploration of complex problems; still, the open challenge for an AI/ML domain is concluded with a short *a priori* knowledge available.^{47} Used in the article, the artificial neural networks (ANNs) exhibit additional advantages; e.g., when trained, a one-step conclusion is realized in the ANN for a given input data, which is of significant meaning for real-time applications, where an iterative scheme of the adaptive procedures can be ineffective; this is the case for a parameter estimation of the time-dependent processes (e.g., MFBM). On the other hand, the challenge is a relatively high expectation of the ANN as regards the amount of data used during training. This triggers fundamental and applied studies in the range of ANN, measurement procedures, and signal processing focused on the development of the methods and tools valid for reliable characterization of complex systems and signals with a limited amount of corrupted data.^{48,49} The problem identified in the paper and the prescription developed for solving it bring contributions both in a range of the studied process (MFBM) and tools (ANN applied to modeling of the fractional dynamics with a limited amount of input data). A recent signal processing literature indicates this research area as very promising, also in the area of fractional dynamics systems.^{50–53}

The other contribution introduced in the paper is a discriminating procedure that can be applied to sort out the trajectories into two classes: classical fractional dynamics systems (the case with a constant Hurst index) and time-dependent fractional dynamics (the case with a time-dependent Hurst exponent). As provided, in the paper, the methodology for intelligent labeling of noisy data with a short history opens new opportunities for signal processing, e.g., toward real-time algorithms for annotation of time series data, where the stated existence of the FBM and MFBM components can indicate time slots with various modes of system operation. This method enables a further practical use of the contribution, e.g., to automatize the monitoring of physiological systems (telemedicine), identification of normal and critical events among the meteorological data, and in industrial applications of IoT (Industry 4.0).

The rest of the paper is organized as follows: In Sec. II, we review the most important properties of the FBM and MFBM. In Sec. III, we introduce the estimation algorithm with the description of the proposed model’s architecture alongside its training methodology and the random dataset construction. Next, in Sec. IV, we check the effectiveness of the proposed technique for the MFBM time-varying output with different Hurst exponent functions. The efficiency of the ANN-based technique is compared with the results of the approach proposed in Ref. 39. In this section, we also propose a step-by-step procedure for discriminating the cases with constant and time-dependent fractional dynamics, characterized by the Hurst index. In Sec. V, we perform the analysis for the experimental data. Section VI summarizes the paper.

## II. FRACTIONAL AND MULTIFRACTIONAL BROWNIAN MOTION

Fractional Brownian motion can be considered an extension of classical Brownian motion.^{14,15} FBM ${XH(t),t\u22650}$ with the Hurst exponent (also called the Hurst index) $H\u2208(0,1)$ is a zero-mean Gaussian process that fulfills the following Langevin equation:^{12,54}

In the above equation, ${\eta H(t),t\u22650}$ is the discrete fractional Gaussian noise (FGN), which is a zero-mean Gaussian process with the following autocovariance function (ACVF):

The parameter $D>0$ is the so-called diffusion coefficient. It is worth highlighting that the FGN is a stationary process. The process could also be defined using fractional integration^{27,37,55} with the Riemann–Liouville fractional derivative,^{56} which is a special case of the Dzherbashian–Nersesian fractional operator.^{57}

The ACVF of FBM is given by

For $H>1/2$, the FGN is a positively correlated process and exhibits long-range dependence (long-memory or persistence).^{12} In this case, the FBM is considered a super-diffusive process. For $H<1/2$, the FGN is a negatively correlated process, and it exhibits short-range dependence (medium dependence or anti-persistence). In this case, FBM is a sub-diffusive process. For $H=1/2$, the FBM reduces to the standard Brownian motion (BM) ${X12(t),t\u22650}$.

The multifractional Brownian motion ${XH(t)(t),t>0}$ is one of the extensions of the classical FBM defined above. The MFBM is a zero-mean Gaussian process with the following ACVF for $t,s\u22650$:^{26,58}

where $D(\u22c5,\u22c5)$ is the function given by

In the above equations, $\sigma >0$ is a constant, $\Gamma (\u22c5)$ is a Gamma function, and $H(\u22c5):[0,\u221e)\u2192[a,b]\u2282(0,1)$ is a Hölder function.

Similar to the classical FBM, the function $H(\u22c5)$ is also called the Hurst exponent. For interesting properties of MFBM, see, for instance, Ref. 37.

In this paper, we consider the following Hurst exponent $H(\u22c5)$ defined for $t\u2208[0,T]$ ($T>0$ is a time horizon)

- Case 1—constant Hurst exponent:(6)$H(t)\u2261H;$
- Case 2—linear Hurst exponent:(7)$H(t)=a1t+a2;$
- Case 3—logistic Hurst exponent:(8)$H(t)=b1+b21+exp\u2061{\u2212b3(tT\u2212b4)};$
- Case 4—periodic Hurst exponent:(9)$H(t)=c1sin\u2061(c2tT)+c3.$

As it was mentioned, case 1 corresponds to the classical FBM. Case 2 corresponds to a process that switches steadily from sub- to super-diffusive regimes (or vice versa). In case 3, we have a similar situation; however, the change takes place more rapidly, which can be used to model a fast transition in the mode of motion.^{59,60} The Hurst exponent considered in case 3 approximates the regime switching Hurst exponent. More precisely, with $b1=H1,b2=H1\u2212H2,b4=T~\u2208(0,T)$, we have

Specifically, the corresponding process approximates a process that switches at time $T~$ from a FBM with Hurst exponent $H=H1$ to a FBM with $H=H2$. A particular case is when $H1<0.5<H2$. This is the situation where the process for times $t<T~$ exhibits sub-diffusive behavior, while for $t>T~$, it is a super-diffusive one. Case 4 corresponds to the situation where the regime changes are gradual and repetitive.

## III. ESTIMATION ALGORITHM

The estimation of the Hurst exponent $H(t)$ for the MFBM based on a single trajectory is a complex problem. The main rationale behind this statement is that we want to estimate ${H(t)}t=1T$ vector of length $T$ given a sample trajectory of length $T$ ${XH(t)(t)}t=1T$. Thus, to estimate the value of the Hurst exponent for a given time $t$, one can rely solely on a considered trajectory at point $t$ and its surroundings. Depending on how locally ergodic the process is, we can use a smaller or bigger window size to estimate $H(t)$. Regarding the limited amount of *a priori* information of the process properties encoded in the input data and the unclear rules governing the selection of the window size for analysis, the problem is solved in the paper with the use of a deep learning methodology based on convolutional deep neural networks.

To enrich the input information about the process, $M=4$ different distance measures have been applied to the raw signal during the preprocessing phase, namely,

${log\u2061(1+|xi\u2212xj|)}ij$,

${log\u2061(1+(xi\u2212xj)2)}ij$,

${log\u2061(1+|xixj|)}ij$,

${log\u2061(1+(xixj)2)}ij$,

where $i,j\u2208{1,2,\u2026,T}$.

Distances (a) and (b) resemble the Euclidean distance, while the last two metrics are in analogy to the covariance. Since the MFBM can be defined through the covariance matrix (like any Gaussian process), the empirical covariance matrix can be used to characterize a given signal. The motivation to use the Euclidean distance comes from the ACVF’s (4) interpretation, i.e., the higher $H(t)$, the higher is the variance of the process and thus raises the probability of increasing the distance between observations. The distance expressed by the function $log\u2061(1+|x|)$ is used to give more significance to values closer to zero, while the distance based on $log\u2061(1+x2)$ is used to emphasize more distant and rapid changes. The $log\u2061(\u22c5)$ function is used to achieve resilience to the outliers and stability of the estimator.

The deep learning model consists of two parts: an encoder and a decoder.^{44} The transformation is made from a three-dimensional ($T\xd7T\xd7M$) input into a one-dimensional vector ($T\xd71=T$) corresponding to the values $(H(1),H(2),\u2026,H(T))$. Both the encoder and decoder are based on Inception-ResNet,^{61} and the whole network resembles U-Net.^{44} The schematic representation of the structure applied to the model is shown in Fig. 1. Details of each block are described in Appendix B.

The selection of Inception-ResNet, as an inspiration for both the encoder and decoder, was motivated by the characteristics of the inception network.^{62} Specifically, combining both more local ($3\xd73$) and wider ($7\xd71$ and $1\xd77$) filters enables to capture a variety of dependencies. Furthermore, this architecture implements residual connections, introduced with the ResNet,^{63} which enables building significantly deeper NNs trained with a backpropagation technique^{64} and improves convergence during the training process.

In the paper, the estimation problem consists of mapping a tensor ($T\xd7T\xd7M$) onto a vector ($T\xd71=T$). As the result, the network (Fig. 1) was also significantly inspired by the U-Net’s architecture.^{44} However, in the classic implementation, the order of the tensor does not change. For this reason, the architecture was modified such that the encoder works on three-dimensional tensors and the decoder on two-dimensional ones. To extract the most important (dominant) feature from each column, a reduction method was applied to the input data where the maximal value over one of the $T$ axes was taken. It is crucial as if the average is taken, and then the information related to the whole signal is used. This stands in contradiction to the locality of the estimates. The network was trained using randomly generated outputs. In this scenario, each time a unique realization of the trajectory is provided to the NN input. The batch consisted of two samples; thus, the input tensor has the following shape: $2\xd7T\xd7T\xd7M$ (as previously stated $M=4$). The $T$ value (corresponding to the number of samples of a trajectory) was randomly selected from $T\u2208{256,384,512}$. If it had not been specified in the other way, the evaluations were based on the largest $T=512$. For a single signal, a random $H(t)$ is selected—one of the functions presented in Eqs. (6)–(9). Moreover, the random Hurst index was added to the learning set, to reconstruct the possible random changes of this parameter, which is especially important during the analysis of experimental data (presented in Sec. V). We include the cases when $H(t)$ is a stochastic process. Here, three processes have been considered, i.e., the reflected Brownian motion (RBM), the reflected Ornstein–Uhlenbeck (ROU) process, and the smoothed telegraph process (ST). The corresponding definitions are presented in Appendix A, and the parameterization for these processes applied during experimentation is specified in Table I.

Case . | Parameters’ sampling method . |
---|---|

Case 1—Eq. (6) | $H\u223cU(0.01,0.99)$ |

Case 2—Eq. (7) | $H0\u223cU(0.01,0.99)$, $H1\u223cmax(min(H0+N(0,0.25),0.99),0.01)$, $a1=H1\u2212H0T$, a_{2} = H_{0} |

Case 3—Eq. (8) | $H0\u223cU(0.01,0.99)$, $H1\u223cmax(min(H0+N(0,0.25),0.99),0.01)$ |

b_{1} = H_{0}, b_{2} = H_{1} − H_{0}, $b3=U(10,250)$, $b4=U(0.1,0.9)$ | |

Case 4—Eq. (9) | $H\xaf\u223cU(0.02,0.98)$, $c1\u223cU(0.01,min(0.99\u2212H\xaf,H\xaf))$, $c2=U(1,200)$, $c3=H\xaf$ |

RBM—Appendix A 1 | $H0\u223cU(0.1,0.9)$, $\sigma \u223cU(0.01,0.99)$ |

ROU—Appendix A 2 | $H0\u223cU(0.1,0.9)$, $\sigma \u223cU(0.01,0.99)$, $\lambda \u223cU(0.1,100)$ |

ST—Appendix A 3 | $H0\u223cU(0.01,0.99)$, $H1\u223cU(0.01,0.99)$, $\gamma \u223cU(0.1,10)$, $\beta \u223cU(0.1,100)$ |

Case . | Parameters’ sampling method . |
---|---|

Case 1—Eq. (6) | $H\u223cU(0.01,0.99)$ |

Case 2—Eq. (7) | $H0\u223cU(0.01,0.99)$, $H1\u223cmax(min(H0+N(0,0.25),0.99),0.01)$, $a1=H1\u2212H0T$, a_{2} = H_{0} |

Case 3—Eq. (8) | $H0\u223cU(0.01,0.99)$, $H1\u223cmax(min(H0+N(0,0.25),0.99),0.01)$ |

b_{1} = H_{0}, b_{2} = H_{1} − H_{0}, $b3=U(10,250)$, $b4=U(0.1,0.9)$ | |

Case 4—Eq. (9) | $H\xaf\u223cU(0.02,0.98)$, $c1\u223cU(0.01,min(0.99\u2212H\xaf,H\xaf))$, $c2=U(1,200)$, $c3=H\xaf$ |

RBM—Appendix A 1 | $H0\u223cU(0.1,0.9)$, $\sigma \u223cU(0.01,0.99)$ |

ROU—Appendix A 2 | $H0\u223cU(0.1,0.9)$, $\sigma \u223cU(0.01,0.99)$, $\lambda \u223cU(0.1,100)$ |

ST—Appendix A 3 | $H0\u223cU(0.01,0.99)$, $H1\u223cU(0.01,0.99)$, $\gamma \u223cU(0.1,10)$, $\beta \u223cU(0.1,100)$ |

Each function has parameters that are selected randomly in the training set using random variables defined in Table I. Then, based on a trajectory, $M$ distance matrices are computed and scaled using the following formula:

Additionally, normalization is applied to the input data:

Parameters $\mu k$ and $\sigma k$ are estimated before the training using $N=1024$ randomly generated samples.

Since the values depend only on themselves, the scaling was applied. This operation enables to control the range of the input regardless of the type of diffusion or $H(t)$. In the study, the first scaling is used to equate the scales of each distance matrix. The input normalization [illustrated as a block of the network in Fig. 13(a)] assures the zero mean and unit variance for the input; this normalization accelerates and stabilizes the convergence of the optimization algorithms.

The training duration was set to 40 epochs, each consisting of 4048 batches. Several methods to monitor and tune the learning process were applied. First, the learning rate scheduler was used, which contains a warm-up phase. The following recurrent function for learning rate profiling was implemented:

On top of that, a learning rate reduction in the plateau (after 3 epochs of stagnation by a factor of 0.3) and early stopping^{65} (after 7 epochs) conditions were added to mitigate jumping around an optimum and to stop the computations when no further improvement is observed, respectively. Both use the validation error to make an informed choice whenever the learning rate needs to be reduced or the training stopped. The Adam optimizer was used to fit the model to the data during training (while the batch size caused by hardware limitations is small, we increased parameter $\beta 1=0.99$).^{66–71} The mean squared error (MSE) function was used as the loss function during optimization.

## IV. COMPUTER VALIDATION OF THE PROPOSED ALGORITHM

In this section, we present the effectiveness of the methodology presented in Sec. III and propose a visual test to discriminate if the real data correspond to the MFBM with a constant or time-dependent Hurst exponent.

The effectiveness of the proposed estimation algorithm is evaluated using $500$ simulated signals of length $T=512$ for cases 1–4 discussed in Sec. II. A time-varying MFBM output is simulated with the use of the algorithm based on the Cholesky decomposition for Gaussian processes.^{72,73} For the validation, we use MFBM with the parameters presented in Table II. The estimation results are presented in Fig. 2. With the light red lines, we present the 95% confidence intervals (calculated based on the obtained estimated Hurst exponent), while with the dark red line, we present the median of $H^(t)$. The blue line corresponds to the theoretical value of the Hurst exponent.

Case . | Parameters . |
---|---|

1(a) | H = 0.2 |

1(b) | H = 0.8 |

2 | $a1=0.8\u22120.21000,a2=0.2$ |

3(a) | b_{1} = 0.2, b_{2} = 0.6, b_{3} = 15, b_{4} = 0.5 |

3(b) | b_{1} = 0.2, b_{2} = 0.6, b_{3} = 250, b_{4} = 0.5 |

4 | c_{1} = 0.3, c_{2} = 10, c_{3} = 0.5 |

Case . | Parameters . |
---|---|

1(a) | H = 0.2 |

1(b) | H = 0.8 |

2 | $a1=0.8\u22120.21000,a2=0.2$ |

3(a) | b_{1} = 0.2, b_{2} = 0.6, b_{3} = 15, b_{4} = 0.5 |

3(b) | b_{1} = 0.2, b_{2} = 0.6, b_{3} = 250, b_{4} = 0.5 |

4 | c_{1} = 0.3, c_{2} = 10, c_{3} = 0.5 |

One can observe that the proposed algorithm estimates $H(t)$ correctly (with wider confidence intervals for small and large times). Both examples—constant functions [Cases 1(a) and 1(b)]—are estimated accurately. The confidence intervals are narrow, and the estimated functions are very close to the theoretical ones. For the linear function (Case 2), one can observe a significant increase of uncertainty at the boundary. This behavior is in line with intuition as there is less evidence that this function continues to be linear outside the boundary. By interpreting the confidence intervals, the model prefers a constant solution at the boundaries. In the following examples, one can observe that the more volatile $H(t)$, the more unstable the estimator $H^(t)$. This can be observed by comparing the twin examples—Cases 3(a) and 3(b). With smoother $H(t)$ [Case 3(a)], the quantile lines of $H^(t)$ are also smoother and less variable compared to Case 3(b), where one can observe sudden jumps and distortions of the estimates. In the situation when $H(t)$ is periodic (represented by Case 4), the estimator underestimates and favors the solution closer to a constant, i.e., the mean value.

We have compared the obtained results using the proposed algorithm with the technique considering the classical method for estimation of the Hurst exponent for MFBM presented in Ref. 39. This method is based on the observation that the convolved trajectory (with a specific wavelet filter) produces a linear relation with the dilatation rate of the filter. In this relation, the slope is proportional to $H(t)$ at a given time $t$, and it can be extracted using basic transformations. However, it cannot be done for all $t\u2208{0,\u2026,T}$ because the filter (wavelet) length amplified by the dilatation limits the number of beginning times $t$ for which the estimate can be evaluated. We have made the comparison using two lengths of signals: $T=256$ and $T=2560$, with the examples proposed in Ref. 39 [with proposed optimal parameters for each of the $H(t)$ functions]. One can observe (see Fig. 3 and the summary statistics in Table III) that for $T=256$, the proposed estimation algorithm outperforms the classical technique. First of all, the proposed method allows one to estimate $H(t)$ for all values of the time points $t\u2208{1,2,\u2026,T}$. This is crucial, especially for small trajectories, because, with such a limitation, we can lose information about a large percentage of a trajectory as the estimate does not exist. The algorithm presented in Ref. 39 (with proposed optimal parameters) is unable to estimate the first 39 values of ${H(t)}t=1T$, effectively producing a vector ${H^(t)}t=40T$.

Case . | T . | Method . | MAE
. | RMSE
. |
---|---|---|---|---|

Linear case | 256 | A | 0.0566 | 0.0732 |

B | 0.0971 | 0.1237 | ||

2560 | A | 0.0380 | 0.0485 | |

B | 0.0560 | 0.0710 | ||

Logistic case | 256 | A | 0.0629 | 0.0800 |

B | 0.1052 | 0.1293 | ||

2560 | A | 0.0340 | 0.0443 | |

B | 0.0549 | 0.0691 |

Case . | T . | Method . | MAE
. | RMSE
. |
---|---|---|---|---|

Linear case | 256 | A | 0.0566 | 0.0732 |

B | 0.0971 | 0.1237 | ||

2560 | A | 0.0380 | 0.0485 | |

B | 0.0560 | 0.0710 | ||

Logistic case | 256 | A | 0.0629 | 0.0800 |

B | 0.1052 | 0.1293 | ||

2560 | A | 0.0340 | 0.0443 | |

B | 0.0549 | 0.0691 |

The comparative studies indicate that the variance of the estimator presented in Ref. 39 is higher than the variance of the introduced estimator. For the case presented in Ref. 39, namely, for $T=2560$, the considered methods coincide with the introduced algorithm having a slightly lower variance. However, it has bias issues on the boundaries. This is mainly caused by the fact that the model was not trained on such long trajectories, and it has learned an estimator $H^(t)$ for which $\u2200t<0H(t)=H(0)\u2227H\u2032(t)=0$ and $\u2200t>TH(t)=H(T)\u2227H\u2032(t)=0$. Despite this, in general, the proposed estimator characterizes the lower error (see Table III with summary statistics); thus, one can conclude that the algorithm presented in this paper is more effective overall, especially, for shorter trajectories.

It is worth highlighting that the algorithm proposed in this paper is superior with respect to the algorithm described in Ref. 39 because it does not need to find optimal hyperparameters. The model is equipped with optimal filters and aggregates the information with an optimal window size, whereas in the method proposed in Ref. 39, it is needed to specify appropriate wavelets, with an optimal aggregation window, where choice has a tremendous effect on the estimate. Choosing a too wide window can lower the variance of the estimator; however, it will increase the bias toward a constant solution. With a too narrow window, the estimator can capture rapid changes in the process’s characteristics; yet, because of the increased variance, it might be random and this could lead to incorrect conclusions. This is difficult, especially because, as it is shown in Ref. 39, the optimal hyperparameters are different for cases with the linear and logistic $H(t)$. This yields that there is a need to know what type of $H(t)$ function the process is governed by before the estimation algorithm starts. Using neural networks, we opt for an optimal solution without the need for guessing the optimal parameter and knowing what type of function $H(t)$ is considered. It is because we train the model with several types of functions and because the convolutional neural networks’ feature is that they have shared filters’ weights. This guarantees for the method to perform well even if a function does not strictly belong to the family of functions the model was trained on. For this reason, the proposed solution will perform better not only for functions originating from training set (see comparison for exemplary functions, with parameterization presented in Table II, shown in Fig. 4 supplemented with summary statistics presented in Table IV) but also for any piecewise-defined function. Furthermore, because we can approximate any function by a piecewise linear function, thus, this method is universal and versatile. We demonstrate it on the example shown in Fig. 5 in which a similar simulation approach is used as in Fig. 3 but with a trajectory length $T=512$ and a more complex $H(t)$ function—the piecewise constant and the sinusoidal function. One can observe that the method proposed in this paper indeed captures the temporal rapid changes of $H(t)$ even though such a function is not in the set of training functions—only its components. Furthermore, from the comparison with the algorithm proposed in Ref. 39, the new method allows one to discover a much more convoluted $H(t)$, whereas the other one smooths out this region causing the temporal periodic behavior to be impossible to educe. This performance difference is even more emphasized in Table V, where the error metrics for Coeurjolly’s method^{39} are almost three times larger.

Case . | Method . | MAE
. | RMSE
. |
---|---|---|---|

1(a) | Proposed method | 0.0277 | 0.0350 |

Coeurjolly’s method | 0.0460 | 0.0578 | |

1(b) | Proposed method | 0.0331 | 0.0420 |

Coeurjolly’s method | 0.0681 | 0.0855 | |

2 | Proposed method | 0.0429 | 0.0551 |

Coeurjolly’s method | 0.0621 | 0.0793 | |

3(a) | Proposed method | 0.0465 | 0.0594 |

Coeurjolly’s method | 0.0684 | 0.0872 | |

3(b) | Proposed method | 0.0436 | 0.0576 |

Coeurjolly’s method | 0.2679 | 0.4248 | |

4 | Proposed method | 0.0693 | 0.0889 |

Coeurjolly’s method | 0.1231 | 0.1475 |

Case . | Method . | MAE
. | RMSE
. |
---|---|---|---|

1(a) | Proposed method | 0.0277 | 0.0350 |

Coeurjolly’s method | 0.0460 | 0.0578 | |

1(b) | Proposed method | 0.0331 | 0.0420 |

Coeurjolly’s method | 0.0681 | 0.0855 | |

2 | Proposed method | 0.0429 | 0.0551 |

Coeurjolly’s method | 0.0621 | 0.0793 | |

3(a) | Proposed method | 0.0465 | 0.0594 |

Coeurjolly’s method | 0.0684 | 0.0872 | |

3(b) | Proposed method | 0.0436 | 0.0576 |

Coeurjolly’s method | 0.2679 | 0.4248 | |

4 | Proposed method | 0.0693 | 0.0889 |

Coeurjolly’s method | 0.1231 | 0.1475 |

### A. Discriminating algorithm between cases with a constant and time-varying Hurst exponent

As it was shown in the above analysis, even if the analyzed data correspond to the MFBM with a constant Hurst exponent (i.e., FBM), the estimated $H(t)$ is not constant. Furthermore, it varies around the true value of the Hurst exponent, which comes from the fact that the calculated $H^(t)$ from the data is only an estimator that, from a mathematical point of view, is a random variable. Hence, there is a need to discriminate the cases when the non-constant $H^(t)$ corresponds to the MFBM with the Hurst exponent being a time-varying deterministic function and FBM [$H(t)$ is constant]. Consequently, in this part, we propose a simple visual test that can be used to discriminate the mentioned cases. The discriminating algorithm proceeds as follows:

For a real signal of length $T$, estimate the Hurst exponent $H^0(t)$, $t=1,2,\u2026,T$ using the proposed estimation algorithm;

Calculate the mean of the obtained values $H^0(t)$ along $t$. This value we denote as $Htest$;

Simulate $M$ FBM time-invariant outputs with the Hurst exponent $Htest$ of the length $T$;

For each simulated signal, estimate the Hurst exponent using the proposed estimation algorithm. As a consequence, we obtain $M$ values of $H^(t)$;

For each $t=1,2,\u2026,T$, calculate the empirical quantiles on the level $\alpha /2$ and $1\u2212\alpha /2$ and construct the empirical confidence interval for the Hurst exponent on the level $1\u2212\alpha $. For a given $t$, we denote it as $[Q\alpha /2(t),Q1\u2212\alpha /2(t)]$;

For each $t=1,2,\u2026,T$ check if $H^0(t)$ falls into the constructed confidence interval. If for minimum $100(1\u2212\alpha )%$ of the $t$ values it is satisfied that $H^0(t)\u2208[Q\alpha /2(t),Q1\u2212\alpha /2(t)]$, then we can suspect that the signal corresponds to the FBM. Otherwise, we cannot assume the constant Hurst exponent and assume the signal corresponds to the MFBM with time-varying $H(t)$.

The discrimination between constant and time-varying Hurst exponents can be performed based on the plot (visual inspection), where we demonstrate the functions $H^0(t)$ and the constructed empirical confidence intervals $[Q\alpha /2(t),Q1\u2212\alpha /2(t)]$ along $t=1,2,\u2026,T$ or on the basis of the percentage of time points for which the estimated Hurst exponent falls into the constructed confidence intervals.

To demonstrate the effectiveness of the proposed discriminating algorithm, first, we analyze the MFBM time-varying output with the linear Hurst exponent corresponding to Case 2, see Eq. (7), and the same parameters used in the simulation study (listed in Table II). Here, we consider the signals of length $T=512$. For the simulated signal, we perform the steps mentioned above, and in Fig. 6(a), we show the estimated Hurst exponent and the constructed confidence interval (and the median) for $\alpha =0.05$. We assume $M=1000$ as the number of FBM time-invariant outputs used to construct the confidence intervals and the median. As one can see, the estimated Hurst exponent does not fall into the constructed confidence interval for most of the $t$ values ($81%$) and the hypothesis of the constant $H(t)$ is rejected.

As a second example, we analyze the FBM time-invariant output with the Hurst index equal to the mean of $H(t)$ used in the linear case described above with the value equal to $H=0.5$. In Fig. 6(b), we demonstrate the visual test described above for $\alpha =0.05$ and $M=1000$. This result indicates that the hypothesis of the constant Hurst exponent cannot be rejected. The percentage of time points for which the estimated $H0(t)$ falls into the constructed confidence interval is equal to $100%$. In this plot, one can also observe that the median obtained from Monte Carlo simulations is not very close to the constant $H(t)\u2261H=0.5$ from which the process was simulated. The reason for that is the time average taken from the series of estimators ${H^(t)}$ resulted in a value only close to the median. This behavior is strictly connected to ergodic theory and the fact that the following limit $limT\u2192\u221e1T\u2211t=1TH^(t)$ does not converge to $H$ as fast as $limM\u2192\u221e1M\u2211i=1MH^i(t)\u223cEH^(t)$ due to the inter-dependencies in ${H^(t)}$ as all of the estimates ${H^(t)}t=1T$ use common trajectory, whereas ${H^i(t)}i=1M$ are evaluated using $M$ independent trajectories.

## V. EXPERIMENTAL DATA ANALYSIS

To evaluate the performance of the proposed algorithm on experimental data, we analyze the trajectories of polystyrene beads of 50 nm in diameter in an agarose hydrogel. The experimental system has been previously described.^{16} A 1.5% agarose gel was prepared from agarose powder by dissolving it in phosphate-buffered saline. The polystyrene beads were first heated to 60 $\xb0$C in 0.5% Tween 20 and introduced into the agarose solution that was maintained at the same temperature and mixed for 15 min. Then, it was transferred to a hot cover-slip where it slowly cooled until it reached room temperature. The beads were imaged in an inverted microscope with a 40$\xd7$ objective at 71 frames/s. Tracking of the individual beads was performed in LABVIEW using a cross-correlation-based algorithm.^{74} A total of 40 trajectories, consisting of 512 data points each, were analyzed here.

It was previously established that the motion of microspheres in agarose hydrogels displays mean squared displacement, power spectral density,^{16} and an empirical anomaly measure^{75} consistent with those of FBM. However, the hypothesis that the Hurst exponent is constant in time has never been rigorously tested. Thus, the MFBM testing algorithm provides a good tool for the validation or rejection of FBM.

Representative trajectories are shown in Fig. 7. To demonstrate how the proposed estimation procedure and the discriminating algorithm described in Sec. IV A can be applied to real data, we demonstrate the estimated $H(t)$ and the constructed confidence intervals calculated for the trajectories of FBM with the Hurst exponent equal to the mean of $H^(t)$ (Fig. 8), similar to what was described in the discriminating algorithm. In the upper panels, we demonstrate the results for the trajectories presented in Fig. 7, for which the estimated Hurst exponent falls into the constructed confidence interval for minimum $100(1\u2212\alpha )%$ time points (see information on the top of the plots), while in the bottom panels, we show the results for selected trajectories for which we have the opposite situation. In the second case, the hypothesis of FBM cannot be accepted, and we assume the MFBM with the time-varying Hurst exponent. Additionally, in Fig. 8, we also demonstrate the median of the estimated Hurst exponent along time points. In the analysis, we take $\alpha =0.05$ and $M=1000$. In general, the number of trajectories for which the hypothesis of FBM cannot be rejected using the introduced discriminating algorithm is equal to $29$, while for $11$ trajectories, the hypothesis of FBM with constant $H$ is rejected.

In Fig. 9, for trajectories presented in Fig. 7, we show $H^(t)$ with the confidence intervals constructed using $1000$ MFBM time-varying outputs with the estimated Hurst index. One can clearly see that the signals simulated with the corresponding models (see confidence bounds) fully reflect the behavior of the estimated Hurst index, such as rapid changes and the time-varying Hurst index for the cases demonstrated in the upper panels of Fig. 9. The cases presented in the upper panels, according to the discriminating algorithm, correspond to the FBM. Indeed, the estimated Hurst indexes as well as the confidence bounds obtained from the simulated signals exhibit small volatility with respect to the cases demonstrated in the bottom panels. Additionally, it should be highlighted that the estimated $H0(t)$ (corresponding to the analyzed trajectories) falls into the constricted confidence intervals. From the above analysis, we can conclude that the tested MFBM models are appropriate for the examined trajectories.

The current results provide a validation for the FBM model in 73% of the analyzed trajectories. However, the classical FBM is rejected in the remaining 27%. The average of Hurst indexes indicates that the trajectories are of the sub-diffusion type, as expected from previous analyses. The deviations from a constant Hurst index are usually observed during specific periods at which the trajectories either become super-diffusive ($H>0.5$) or are seen to exhibit a more marked sub-diffusive behavior. In the former cases, a plausible explanation is that small transient water flows can cause the beads to exhibit directed motion such that the motion appears super-diffusive during these short times. On the other hand, the case of a reduced Hurst index can be caused by imperfections in the hydrogel leading to material heterogeneities. In particular, as a particle traverses a denser region within the medium, it would experience more negative increment autocorrelations and it would transiently appear as FBM with smaller $H$.

While the majority of the trajectories are found to be of the FBM type, validating previous assumptions, the use of the algorithm for the identification of MFBM in individual trajectories provides very useful information along two different lines of data analysis. First, the method can be used to discriminate trajectories that do not correspond to the simple FBM with constant $H$. Furthermore, the evaluation of the deviations from constant $H$ provides a tool for the characterization and understanding of these anomalous events. Second, given that a viscoelastic material with homogeneous physical properties should yield trajectories with constant $H$, the detection of MFBM can be used to evaluate the properties and quality of the material itself. This is, in fact, one of the main goals of microrheology experiments, where the motions of embedded particles are analyzed in order to measure the mechanical properties of the medium on microscopic scales.^{76}

## VI. SUMMARY

The contributions provided in the paper can be listed in several highlights:

The problem of identification and characterization of the complex fractional dynamics is discussed and solved for limited input information and the exemplary model of MFBM [with time-varying $H(t)$ coefficient, adequate to the temporal mechanisms of the fractional system dynamics].

The AI-based scheme has been designed for constant and time-dependent $H(t)$ estimation and for the sequential binary annotation of the measured signal.

The $MAE$ and the $RMSE$ estimation errors, measured during the computational experiments for linear, piecewise, and sinusoidal changes of $H(t)$, did not exceed 6% and 8%, respectively, which is half that of the reference method.

The methodological prescription outlined qualitatively and quantitatively in the paper can be mapped to complex systems of any nature where the observed output exhibits temporal fractional dynamics.

The designed AI-based algorithm can inspire and contribute to other real-time signal processing schemes operating with short time series in the presence of noise. This work paves the way for developments in various services, including, e.g., telemedicine, IIoT/Industry 4.0, smart city, etc.

## ACKNOWLEDGMENTS

The work of A.W. was supported by the National Center of Science under Opus Grant No. 2020/37/B/HS4/00120 “Market risk model identification and validation using novel statistical, probabilistic, and machine learning tools.” D.K. acknowledges the support from the National Science Foundation under Grant No. 2102832.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Dawid Szarek:** Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (lead); Validation (equal); Visualization (lead); Writing – original draft (equal); Writing – review & editing (equal). **Ireneusz Jabłoński:** Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Project administration (lead); Supervision (lead); Writing – original draft (equal); Writing – review & editing (equal). **Diego Krapf:** Conceptualization (equal); Formal analysis (equal); Investigation (equal); Resources (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). **Agnieszka Wyłomańska:** Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Project administration (lead); Supervision (lead); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX A: STOCHASTIC PROCESSES CORRESPONDING TO THE HURST EXPONENT USED IN THE LEARNING SET

#### 1. Reflected Brownian motion

Brownian motion ${HBM(t),t\u22650}$ is a stochastic process defined by the following Langevin equation:^{77}

where ${\eta 1/2(t),t\u22650}$ is the Gaussian white noise process (i.e., process of uncorrelated random variables with zero-mean and constant variance—here, we assume that it is equal to $1$) and $\sigma >0$. The ACVF for BM is given by

The reflected Brownian motion is a way of keeping the Brownian motion within a given domain. In our case, the domain is given by a $(0,1)$ interval, and the reflection is performed at any given time when the process exceeds it. In the general case, it can be any given domain.

The process we use for $H(t)$ that differs from the definition (A1) is that the starting point $HBM(0)=H0$.

#### 2. Reflected Ornstein–Uhlenbeck process

The classical Ornstein–Uhlenbeck process ${HOU(t),t\u22650}$ is a stationary solution of the Langevin equation of the following form:^{78}

where ${\eta 1/2(t),t\u22650}$ is the Gaussian white noise and $\mu \u2208R$, $\lambda >0$ and $\sigma >0$ are the parameters of the process. On the other hand, the OU process can be obtained from a Brownian motion by the so-called Lamperti transformation.^{79}

When $HOU(0)$ is constant and $\mu =0$, the ACVF of the considered process is given by

The reflected Ornstein–Uhlenbeck process is defined similarly as the reflected Brownian motion defined in (1). Specifically, at every given time $t$, when a trajectory exceeds a given threshold, it is reflected.

The parameterization for trajectories’ simulation is that $\mu =H0$ and $HOU(0)=H0$.

#### 3. Smooth telegraph process

The telegraph process (asymmetric) ${HT(t),t\u22650}$ is defined as follows:^{80}

where ${N(t),t\u22650}$ is the homogeneous Poisson process^{81} with the intensity $\gamma >0$ and $h0,h1$ are the positive real numbers such that $h0<h1<1\u2212h0$. The asymmetric telegraph process is a two-stage process taking the values $H0=h1\u2212h0$ and $H1=h1+h0$ with probabilities $e\u2212\gamma tsinh\u2061(\gamma t)$ and $e\u2212\gamma tcosh\u2061(\gamma t)$, respectively. The smoothed (filtered) asymmetric telegraph process ${HST(t),t\u22650}$ is defined through the following stochastic differential equation in the following form:^{80}

where $\beta >0$ is a parameter of the filter and ${HT(t),t\u22650}$ is an asymmetric telegraph process defined in Eq. (A5).

### APPENDIX B: DETAILS OF NEURAL NETWORK’S ARCHITECTURE

Building blocks of the architecture presented in Fig. 1 are described in detail in the following figures:

Fig. 13(a)—Stem,

Fig. 10(a)—Block A,

Fig. 10(b)—Reduction Block A,

Fig. 11(a)—Block B,

Fig. 11(b)—Reduction Block B,

Fig. 12(a)—Block C,

Fig. 12(b)—Transpose Block C,

Fig. 13(c)—Max-Reduce Block,

Fig. 11(d)—Transpose Reduction Block B,

Fig. 11(c)—Transpose Block B,

Fig. 10(d)—Transpose Reduction Block A,

Fig. 10(c)—Transpose Block A,

Fig. 13(b)—Transpose Stem.

## REFERENCES

*2000 IEEE International Conference on Acoustics, Speech, and Signal Processing. Proceedings (Cat. No. 00CH37100)*(IEEE, 2000), Vol. 6, pp. 3810–3813.

*2020 IEEE 15th International Conference on Industrial and Information Systems (ICIIS)*(IEEE, 2020), pp. 174–179.

*International Conference on Medical Image Computing and Computer-Assisted Intervention*(Springer, Cham, 2015), pp. 234–241.

*et al.*, “

*Thirty-First AAAI Conference on Artificial Intelligence*(AAAI, 2017).

*Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition*(IEEE, 2015), pp. 1–9.

*Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition*(IEEE, 2016), pp. 770–778.

*Proceedings of the IEEE International Conference on Computer Vision*(IEEE, 2015), pp. 2758–2766.

*Proceedings of the IEEE International Conference on Computer Vision*(IEEE, 2017), pp. 2794–2802.

*COMPSTAT*(Springer, 1998), pp. 233–238.