Estimating correlation dimension of high-dimensional signals - quick algorithm

Distinguishing between random stochastic processes and that being highly complex, however regulated by determined differential equations and possessing some kind of hidden regulation is one of the mathematical problems being analyzed in different fields of science. The problem is especially complicated in the case of signals that are suspected to undergo regulation, however, the variability of the process looks to be random. The estimation of the complexity becomes difficult in the case of complicated systems in which many factors interact with each other in non-linear way. One of the proposed methods describing the complexity of the system is the estimation of the correlation dimension d of the signal. The first algorithm was proposed by Grassberger and Proccacia,¹ however, it has been shown to be ineffective for HDS. The algorithm that possesses the highest power for estimating this parameter is the Takens-Ellner algorithm^4,3,2 which was modified by Michalak^5,6 for the purpose of high-dimensional analysis. These papers showed the possibility to estimate the correlation dimension up to about d=8 and probably more.

Estimating correlation dimension d consists of a few steps which are the separate problems from the methodological and algorithmic point of view. The methodological aspects focus on the correct and precise d estimation while the algorithmic ones focus on speeding up the calculations and reducing the memory space used while calculation. Both criteria are as a rule in opposition to each other because the small improvement of the precision forces as a rule huge increase in calculation time and memory space. Keeping both high precision and acceptable calculation time is often a “balance on the line.” Thus, in order to properly use the presented method, all the parameters determining speed and time must be perfectly understood.

Current paper shows brief recall of the algorithm presented in Refs. 5 and 6 and illustrates how to speed up the algorithm without significantly affecting its accuracy.

Current paper focuses on the following aspects of the algorithm:

1. Presenting new ideas improving and accelerating the algorithm.

2. Improving the precision of the elementary Takens-Ellner formula.

3. Defining and reduction of the problem of δ-vector reduction.

4. Presenting the problem of strongly non-Gaussian signals and showing some aspects connected with the preliminary normalization of the signal.

5. Presenting the example of the signal possessing 2 plateaus in the relation $dw$=fn(W).

The  Appendix presents the use of the free available Matlab toolkit with the implementation of the algorithm. This chapter is dedicated for scientists who want use all the implemented options of the algorithm. The HDS-toolkit is available on author’s website “www.drmichalak.pl/chaos/eng/”.

The current section recalls the main points of this algorithm.

Briefly, the algorithm consists of the following steps:

• Choosing the set of window widths W for embedding procedure.

• Estimating $dw$ for every W.

• Fitting some functions (e.g. 4-th degree polynomials) to $dw$=fn(W) and looking for lowest slope point in the fitted function to be the base for estimation of the correlation dimension of the signal.

• Comparing the results with that obtained for surrogates.

Let us analyze and discuss the main problems connected with subsequent steps:

ad I. Choosing a set of 30-50 values of window width W ranging approximately from 0.1τ to 10τ (τ is the time at which the autocorrelation time of the signal drops to 1/e). W denotes the length of the extract of the signal creating one m-dimensional point in the phase space.⁷ Let us denote $dw$ to be the correlation dimension estimated using given W. Strong attention must be paid for signals possessing pink noise-like spectrum, because such signals may possess a very large autocorrelation time. Estimating correlation dimension may be ambiguous in this case. In the case of low-dimensional signals or signals possessing different features in different time scales, it may be necessary to perform the analysis for longer W-values.

ad II. Estimating $dw$ for every W. This step is the main point of the algorithm which must be optimized according to 2 opposing criteria: precision and computation time. It can be divided into a few separate points:

• A.

  Choosing embedding dimension m. The minimum required m is m*=2$dw$+1. It should be, however, about 3-4x larger than the expected estimated $dw$⁠. $dw$ can be expected to be about 5-20% higher than that estimated for the former, a bit smaller W.

• B.

  Embedding the signal using cubic interpolation, choosing the lag L = W/(m-1). Use the cubic interpolation for non-integer L. The choice of J (distance between the beginnings of signal extracts building subsequent embedded points) is not crucial for the precision, however, it can be crucial for memory use because it determines the size of the matrix A[n, m] with embedded points (n ≈ N/J).

• C.

  Calculating the distances between embedded points (building the vector δ) and sorting it in ascending order. Repeating the Random Joining Procedure (RJP) until the assumed number of distances D_a between embedded points is reached. This is the main time-consuming point of the algorithm. The time T of distances’ calculation depends mainly on the embedding dimension m used (T_m∼m) and on the correlation dimension of the signal (T_d∼n^d) (due to the relation log(P_max)∼d presented in Ref. 6). The time for δ-sorting procedure depends on the length of vector δ: T_s∼n · log(n).

• D.

  Building the Correlation Integral function (C’) and looking for left-side plateau in this function.

• E.

  Estimating $dw$ using the classic Takens-Ellner formula 1 for the set of P_max values corresponding to the left-sided plateau. The mean over this range can be treated as the final estimator of $dw$⁠. See Figure 3. (I_min/max = round(P_min/max · length(δ)):

• F.

  Repeating the iteration using higher m if the estimated $dw$ is higher than (m-1)/2 in order to fulfill the rule m > 2$dw$+1. (In practice, I propose to use the saturation reserve $dw$ < m/3 in order to estimate the correlation dimension of surrogates using the same m which makes it possible to observe the increase in $dw$ using exactly the same estimation parameters.)

ad III. This step does not consume a significant amount of time. The question related to this point concerns the choice of the type of fitting function used to find the Lowest Slope Point in $dw$=fn(W) relation which determines the final estimation of d_opt. The 4-th degree polynomial fitting (see Figure 4) was used in my former paper, however, other fitting functions are possible, as well, and they may give some additional benefits. The use of higher degree polynomials and custom-defined functions is prepared in my HDS-toolkit.

The second question concerns the problem, which P_max should be used to estimate individual $dw$⁠. In my former papers, I proposed to estimate the 2-dimensional grid of $dw×Pmax$ values in order to fit first the set of polynomials to $dw$=fn(W) for individual values P_max. In the current paper, I propose to estimate the final $dw$’s for individual W’s and to perform one final fitting to $dw$=fn(W) relation which is more correct with respect to the concept of the necessity to consider the deflection error while estimating $dw$⁠.

ad IV. Comparing the results with that obtained for the surrogates.^8,9,10 The increase in d after signal shuffling points to some kind of the hidden chaotic regulation. The increase in d is especially expected in the range of plateau/“lower slope area”(LSA) of the relation d=fn(W). In the case of ambiguous results concerning the determination of the plateau in $dw$=fn(W), the point of maximum difference between the surrogates and original signal can be a hint for the determination of the proper plateau region and estimating the final d_opt.

A little problem occurs due to the fact that the correlation dimension of the surrogate is infinite in theory which can lead to the looped repeating of the calculations using increasing m while the estimated d is larger than (m-1)/2. In practice, the calculations are recommended to be performed using the same parameters (W, P_min, P_max, D_a, m, J) as for the original signal in order to obtain direct comparison of $dw$ calculated for the original signal and its surrogates. The embedding dimension m used for the original signal must be only taken with some additional reserve to observe $dw$ increase for surrogates.

The algorithm presented in Ref. 6 does not estimate the optimum $Pmaxo$ and δ_max and estimates $dw$ for a number of different P_max’s up to the constant large value P_max ≈ 0.1 for all $dw$’s being estimated. The full relation d=fn(W, P_max) can be built in this case, to be a grid fulfilling all the 2-dimensional area W x P_max. In the case of large W’s for which large $dw$’s are estimated, the plateau range P_min − P_max ends at very small P_max’s of about 10⁻⁵ − 10⁻³. It means that the large majority of distances, more than 99.9%, is calculated unnecessary. Thus, the simple way to reduce the computation time is to remove these distances from the δ-sorting procedure. The possibility to estimate δ_max to be the distance representing the approximate end of plateau range, makes it possible not to include the longer distances to the vector δ. The final vector δ has the length of about 10³ − 10⁴ and not 10⁶ − 10⁸ in this case.

One of the important conclusions presented in Ref. 6 is the possibility to estimate the approximate end of plateau in $dw$=fn(P_max) relation based on the assumed deflection error e in the histogram convolution plot. The approximate end of plateau is represented by 2 values: δ_max and optimum P_max (⁠$Pmaxo$⁠). Let us recall that the distance δ_max depends only on the used embedding dimension m.⁶ On the other hand, $Pmaxo$ depends on $dw$ to be estimated.

δ_max can be estimated for the assumed deflection error e from the correlation integral $Cm′$ relation being built from the distribution (histogram) of the differences between signal samples.⁶ (First, the histogram H₁ of samples’ differences is calculated using equal bin widths. H_a+b is calculated by a special kind of square root convolution from H_a and H_b staring from a=1 and b=1 and repeating it until the used embedding dimension m is reached. Next, cumulative sum of H_m (⁠$Hmc$⁠) is calculated. $Cm′$ is the local slope of $Hmc$ in the logarithmic scale.) Estimating $Pmaxo$ and δ_max makes it possible to reject longer distances giving significant acceleration of the algorithm.

As presented in Ref. 6 (Section VI), when analyzing the distances between random points using given embedding dimension m, the slope of the cumulative sum $Cm′$ of the distribution of distances in log-log scale deflects gradually for increasing P_maxs from m to 0. Higher value P_max denotes higher δ_max used and bigger ratio of distances included to δ which strongly accelerates the algorithm. The calculations are performed, however, with a bit lower accuracy, due to a bit higher deflection error e corresponding to higher $Pmaxo$ and δ_max. The exemplary value e=5% denotes the possibility to down-estimate $dw$ by 5%. It is presented that δ_max depends on e and m and can be estimated for assumed e before calculations. Thus, all the distances longer than δ_max can be immediately rejected from the vector δ and not be included to the δ-sorting procedure. Additionally, the 2-step rejection is implemented in order to additionally reduce the computation time. First, the distances based on the first and last dimension only of the m-dimensional points are calculated and rejected if the distance is longer than δ_max. This step rejects large percent of long distances. Next, the distances are calculated based on all dimensions and the remaining longer distances are rejected. This procedure is especially effective in the case of large embedding dimensions used for large W’s and consuming majority of calculation time.

In practice, δ² can be kept in memory to reduce the calculation time for square root operation. d estimated using δ² is just exactly 2x bigger than the real one.

It is worth mentioning that the nonlinear relation δ_max=fn(m) presented in Figure 9 in Ref. 6 is (nearly) linear for normally distributed signals for $δmax2$ (see Figure 1) which makes estimation of δ_max from m and e very easy. In the case of preliminary normalization of the signal, the ready to use values from Figure 1 can be used. (They are stored in deltamaxM_PmaxM_normal.mat file in the HDS-toolkit). In another case, such relations must be built individually for every signal for every new m used in subsequent W-steps. The normalization of the signal accelerates in this way the algorithm, because the convolution procedure becomes unnecessary.

$Pmaxo$ can be estimated for the given e, as well. $Pmaxo$ decreases strongly with the correlation dimension of the signal. In case of large $dw$⁠, the small changes in assumed e cause the relative large changes in P_max to be accepted. E.g. in the case of $dw$=16 (as for e.g. W ≈ 10τ), $Pmaxo$ corresponding to e=3% is equal to 10⁻⁹, that corresponding to e=5% is equal to about 6 · 10⁻⁸ and for e=10% P_max is 5 · 10⁻⁶ (see Table I and Figures 8, 9 in Ref. 6). This example shows that increase in the assumed deflection error from 3% to 5% increases 60x the number of distances D to be included to calculation (and correspondingly from 3% to 10% - 6000x). In other words, in case of e=3%, about 10¹¹ distances must be calculated in order to find only 100 being accepted for $dw$ estimation.

The more general estimation of the gain of computation time can be made based on the Figure 8 and equations from Table I in Ref. 6. Let us assume that the computation time is proportional to the inverse of $Pmaxo$⁠. The gain G in computation time when increasing the deflection error from e₁ to e₂ (e.g. from e₁=3% to e₂=5%) can be approximated as:

where a₁, a₂ are the directional coefficients of subsequent straight lines presented in Figure 8 and Table I in Ref. 6. Parameters b, as close to 0, can be omitted in rough G estimation.

It should be pointed, however, that the estimated $Pmaxo$ can differ from the real P_max which reflects the ratio of the number of distances shorter than δ_max and the total number of distances analyzed.

Let us come back to the algorithm. The first step of $Pmaxo$ and δ_max estimation consists in the building of the histogram convolution plots up to the embedding dimension m used for the given W iteration. Next, the acceptable deflection error must be defined. It is proposed to use e=3% for the initial W’s and gradual increase in e for increasing $dw$’s being estimated for subsequent W’s if the calculation time becomes too long. The small increase in the acceptable e gives in this case the significant reduction of the computation time.

In the case of very large $dw$ being estimated as a rule for large W, the use of small deflection error e=3% makes it very difficult to find distances that are shorter than estimated δ_max. Thus, in order to reduce computation time, it is proposed to increase e for subsequent W’s to keep rather fixed computation time and only slightly reduce the precision. Estimation of $dw$ becomes possible in acceptable, predictable time, in this way. The assumed deflection error can next be used to correct of the estimated d=fn(P_max) relation: d_corrected = d_estim · (1 + e/100). The parameter MaxIterationTime=300s was used in the example presented in Figure 2b which forced the increase in e by 1% in the next W-iterations if this time was exceeded.

(In the example from Figure 2 the signal being the sum of 3 Lorenz signals possessing the time scales ratio 1:$2$⁠:$3$⁠, thus possessing the correlation dimension of abut d ≈ 6, was analyzed. The Lorenz system is described by means of 3 equations: dx/dt=10z(y - x), dy/dt=-xy+28x - y, dz/dt=xy-(8/3)z. The variable x was taken to the analysis. The individual Lorenz attractor possesses the correlation dimension that is equal to about d=2.017.)

The example from Figure 2 shows the old (full, Fig. 2a) and new (quick, Fig. 2b) approach. The relation presented in Figure 2a and in my former paper⁵ have been built for all the values P_max ranging from nearly 0 to nearly 1. This approach was connected with the calculation of a large number of distances stored in δ which forced the huge computation time. However, it made it possible to fit the 4-degree polynomials for every P_max used. Now, the new proposed approach is connected with rejection of all distances larger than previously estimated δ_max as not forming the plateau in $dw$=fn(P_max). This approach does not make it possible to build full relation d=fn(P_max, W) (see Figure 2b) and requires building one relation $dw$=fn(W) which is the base for the final polynomial fitting procedure (see Figure 4).

Two regions can be shown in Figure 2b. That for small W’s was estimated using e=3%. The estimated $Pmaxo$ decreased while $dw$ and the computation time increased. The used P_min value is higher for small W’s because a small total number of distances was enough to perform $dw$ estimation. I_min=1 corresponds to relative large P_min. In the second region, $Pmaxo$ does not change significantly because e increases slowly in subsequent W iterations in order not to elongate the computation time.

The algorithm seems to be easy, however, in order to properly estimate d the attention must also be paid to the following problems:

1. The proper preparation of the signal (especially: sampling frequency and noise reduction). The sampling frequency should be large enough to obtain the autocorrelation time of the signal τ higher than 20 samples. If τ is lower then the cubic interpolation error can dangerously increase.

2. The choice of the set of window widths W to be used to create the relation $dw$=fn(W).

   The larger number of W’s being used increases the precision of the relation $dw$=fn(W). About 30-50 values W are proposed be chosen from the range 0.1-10τ. A special attention should be paid to the range 2-4τ representing the range of the expected d_opt. Large values W corresponding to about 8-10τ are connected with a huge part of computation time, because the estimated $dw$ is high in this case (even up to $dw$=20-30) forcing the use of huge embedding dimension m=50-80. So high m forces the huge number of distances to be calculated, because, as presented earlier, the optimum P_max is very small. It happens that after a few minutes of calculations no distances shorter than assumed δ_max are found making the estimation impossible. The question must be answered in this case: (a) Follow calculations? (b) Increase e and δ_max (reduce precision but increase the chance to find some distances to estimate $dw$? (c) Abandon calculation for given $dw$? The automatic increase in e and δ_max being implemented in HDS-toolkit seems to be the best solution of this problem.

3. The choice of the set of values P_min and P_max defining the range of vector δ used to estimate $dw$⁠.

   Vector δ stores the distances between embedded points and represents the ‘inverse of Correlation Integral’, as presented in Ref. 6. The optimum P_min value depends especially on the noise contamination/reduction. The value $Pmaxo$ depends on the assumed precision e and on the value $dw$ to be estimated. Thus, from the formal point of view, the choice of $Pmaxo$ is possible AFTER $dw$ is measured. The following approach is proposed to omit this loop:

   • estimate preliminary $Pmaxo$ using $dw$ from the former W iteration.

   • estimate d_init using this $Pmaxo$⁠.

   • again estimate preliminary $Pmaxo$ and δ_max using d_init.

4. The choice of P_min is mainly determined by noise contamination and sampling precision. The simplest way is to take I_min=1 and not to reject any shortest distances from analysis. If P_max is large enough, this choice does not introduce significant error to final $dw$ estimation. However, it is recommended to reject the left range of $dw$=fn(P_max) relation in which strong oscillations are visible (see Figure 3). I propose to use as P_min the smallest value P_max for which estimated $dw$ differs from plateau level by less than 30%.

5. The second question concerns the method of estimation of plateau level in $dw$=fn(P_max) (see Figure 3). The problem can arise if the large lacunarities occur in this relation or there is a lack of horizonal plateau. Such lack of horizontal plateau may be a feature of stochastic signals, a consequence of non-reduced noise or just a feature of the given signal. I propose to estimate $dw$ as a mean over d’s estimated for P_max values ranging between 1.2·P_min and $Pmaxo$ using logarithmic scale. However, if the range $Pmin−Pmaxo$ is very narrow, the only $Pmaxo$ value can be taken, as well.

6. The choice of the method for estimation the final d (d_opt) from the relation $dw$=fn(P_max, W) (full method) or $dw$=fn(W) (quick method).

   In the former papers^6,5 the 4-th degree polynomial fitting was proposed to be applied to the relations $dw$=fn(W, P_max) estimated for individual combinations P_min/P_max. The lowest slope point of this polynomial (W_opt, d_opt) corresponds approximately to the linear scaling region being observed for low-dimensional signals and was proposed to be the point for reading the value of the correlation dimension d.

   The 4-th degree polynomial is characterized by one minimum only in its 3-degree derivative. Thus, it can give the unambiguous result. The shape of the relation $dw$=fn(W) can be, however, in some cases so untypical that this function may give unprecise fitting. The problem of analysis of the relation $dw$=fn(P_max, W) is, in general, open for individual solutions depending on the individual shape of this relation. Currently, some additional fitting functions have been added to the algorithm. The file “optimal-fitting.m” makes it possible to add easy other functions to solve individual problems.

7. The introduction of e-correction to the relation $dw$=fn(P_max).

   The relation $dw$=fn(P_max) is charged with the error connected with the use of given P_max. Thus, this relation can be increased by e% (see Figure 3) in order to estimate more precisely $dw$’s and next d_opt. This error, is however, as being estimated using random data series, only a rough estimator of the real error reflecting the distribution of points of the signal. Because of that, I leave decision concerning the use of this correction for the HDS-toolkit users. The analytical approach to this problem is recommended, as well.

8. Determining the cut-off distance δ_cut while calculating the distances δ.

   This parameter is similar to δ_max but denotes the distance that is a real cut-off limit, while δ_max denotes the limit for $dw$ estimation. The general question arises in this point: Does it make a sense to use δ_cut being larger than δ_max? It should be pointed that δ_max and $Pmaxo$ estimated statistically using assumed deflection error e may not correspond to real optimum δ_max and $Pmaxo$⁠. Keeping a bit more distances in memory than estimated δ_max may avoid repeating calculations when the estimated δ_max appeared to be a bit to small or the number of distances D_a was not reached in the assumed time which requires repeating iteration using higher e. The margin distances higher than δ_max can be introduced to the sorting procedure and further calculations if necessary saving the time for the iteration repetition.

The current section focuses on the precision of the elementary Takens-Ellner formula (Eq. 1) estimating d using a set of distances from δ ranging from I_min to I_max. As originally presented, the sorted vector δ, when presented in relation to its index, is linear in log-log scale and the directional coefficient of LSR is equal to a=1/d. Therefore, let us consider the artificial vector δ=c(1, 2, 3,…)^1/8 that corresponds to the ideally distributed distances between points of phase space. The estimated d_estim should be precisely equal to d_real=8 for every range of this vector (for every combination of I_min and I_max). The use of the classic formula (Eq. 1): d_estim=D/(Y +(I_min·y(0)) does not return the value d_estim=8, however, but it provides an underestimated value of d. It means that this formula is imprecise in estimating d. This observation was pointed out in Ref. 3 and treated as a fault of the algorithm. Following this observation the error Δd =d_real − d_estim, presented as function of I_min and I_max was calculated. The error values are presented in Figure 5a.

It can be observed that this error is rather small for large I_min and I_max and becomes significant for small I_min/I_max. The error has been observed to be proportional to d_real and does not depend on the scale factor c. It has been found that the relation from Figure 5a creates the oblique plane in log-log-log scale. Thus, after the plane equation was fitted to the obtained error points, the approximate equation for d-error has been empirically obtained. It is presented in Eq. 3:

The detailed analysis of the TE formula showed that the presented error is mainly connected with the imprecise reference value “I_min·y(0)” which determines the initial point of LSR (see Eq. 1). Because of that, it is proposed to change this part of TE formula to:

This change caused the great reduction of the generated error Δd. The repeated error analysis showed that the relation possessed again the shape of the plane in log-log-log scaling. The significanlty lower error, after plane equation fitting, could be now estimated using formula:

For example, if I_min=1 and I_max=10 then the error is equal to Δd=0.958 if the original “I_min·y(0)” is used and the error is equal to Δd=0.064 for the corrected formula 4.

In the last step, Equation 5 must be converted to calculate d_real from d_estim. It is presented in Equation 6:

Now, let us consider another kind of methodological error. First, let us define, as in the former section, the ideally distibuted artificial vector with distances δ=c · i^a (a=1/d and i=1, 2, 3…). This vector can be presented in the form: log δ = a · log(i) + log c, as well, showing the linear scaling in log-log scale. Now, let us take every r-th value of this vector (r - reduction factor). It can be observed that the scaling proportional to a = 1/d can only be obtained if the reduction starts from the index q=r: δ_r = δ[r, 2r, 3r, …]. The formula defining the reduced vector δ_r = c · (r · i)^a can be converted to: log δ_r = a log i + (a log r + log c), where b = a log r + log c denotes the constant value of the straight line equation possessing the directional coefficient a. Let us observe that such reduction of vector δ does not change the estimated value d.

Next, let us analyze another kind of δ-reduction. Let us take every rth value staring from q < r. Now, the vector δ_r is: δ_r = δ[q, r+q, 2r+q, 3r+q,…]. It can be also presented as: δ_r = c · (r · i + q)^a or log δ_r = a log(r · i + q) + log c, (i=0, 1, 2, 3,…). Now, the precise linear scaling is disturbed, because the logarithm of the sum cannot be reduced, and the estimated d is lower.

Figure 6 shows the exemplary relation between d_estim and initial index q varying from 1 to r for r=5, 7, 9,…, 19 and for the exemplary parameters I_min=2 and I_max=50. d was estimated using the modified TE formula presented in Eq. 4 and with the correction from Eq. 6. As expected, the correct d=8 is only observed for q=r and d decreases nearly linearly for decreasing q.

Now, let us observe that the random choice of distances from the ideally distributed vector δ, which happens while practical d estimation, corresponds rather to the regular reduction r starting from the index q=r/2+0.5. E.g. if i=1000 and r=5 than the reduced vector δ_r should start from q=3 (δ_r=δ[3, 8, 13,…, 998] in order to be centrally distributed in δ, as in the case of random choice. The estimated d is lowered in this case and it is a source of methodological error. The difference between d_real=8 and d estimated for q=r/2+0.5 could be treated as an initial estimator of this error. The comparison between d_q=r/2+0.5 and mean d, obtained by random choice of distances from δ calculated as geometric mean for j=10000 repetitions, does not, however, confirm this hypothesis. The values of d_estim obtained by random choice are higher and correspond approximately to the choice of q ≈ 0.7r rather than to q=r/2+0.5 (see Figure 6). Similar results were obtained for other combinations of I_min, I_max and r. This phenomenon has not been explained at the moment and it can be probably explained using analytical methods. Currently, an approximate empirical correction is proposed. As in previous section, it has been observed that the error △d=d_real − d_q=0.7r depends nearly linearly on I_min and I_max in a 3 × log scale. The plane equation fitted to △d=fn(I_min, I_max) in this 3 × log scale gives the approximate equation presented in Eq. 7:

After conversion, d_real can be estimated using the formula in Eq. 8:

The Equations 7 and 8 are precise for r > 5 but small inaccuracies can occur, especially for r < 3. The problem that remains unsolved is the estimation of the value r which is proper for the given vector δ estimated for the given signal, because the distances in δ can be treated as a subset of all the possible distances capable of creating the ideally distributed δ. If r could be estimated then Equation 8 could be slightly modified thereby increasing precision through only a little.

To sum up, the final correction of d requires two corrections - first, that from Eq. 6 and then, that from Eq. 8. Both these corrections are implemented in author’s HDS-toolkit.

The analyzes presented in former paper showed the results for the Gaussian distributed random data series. Now, as an opposite, let us present a short analysis of some strongly non-Gaussian distributed random data series to show significant differences between them. The analyzed series was created from the uniformly distributed series x ∈⟨−1, 1⟩ converted using equation x_ng = sign(x) · (1 − x²) + 1.

Figure 7a shows the distribution of x_ng, distribution of distances (differences) H₁ and slope $C1′$ being the derivative of the cumulative sum of H₁ (see Eq. 3 in Ref. 6). It can be observed that the distribution of distances H₁ does not possess any left-sided plateau. It causes the slope $C1′$ not to be equal to 1 for small distances but to be reduced to about 0.7-0.8. The range of distances, for which the slope is equal to 1 is reduced practically to the point (0, 1) making the estimation of d practically impossible. Figure 7b shows the self-convolutions of H₁ (details of the algorithm are presented in Ref. 6, page 5). The $Cm′$ relations estimated for higher m’s show the strong tendency to down-estimate d.

Basing on the results presented, the conclusion can be drawn that the width of LSR depends mainly on the width of plateau in the H₁. The plateau occurs in the same range of distances for consecutive embedding dimensions. The problem occurs, however, what could be the distribution of the signal that could give the maximal width of left-sided plateau in H₁? The comparison of the normal distribution with a set of other distributions shows that the normal distribution is favorable from this point of view suggesting every non-Gaussian distributed analyzed signals to transform preliminary to the Gaussian ones in order to increase P_max used and reduce in this way the computation time.

The example from the former section discovers the problem of the preliminary normalisation of the signal being analyzed. Normalization is meant as the nonlinear monotonic transformation that causes the distribution of the signal to be normal. Should we normalize the signal before estimating its correlation dimension?

There is a significant benefit of normalisation. Normally distributed signal possesses the broadest plateau in H₁ relation of its distribution. The relative small improvement for low embedding dimension becomes significant for high and very high embedding dimensions used which contributes to significant reduction of computation time. Let us see Figure 8ab which shows the effect of normalisation of the signal being the sum of 3 Lorenz signals possessing the time scale ratio 1:$2$⁠:$3$⁠. The distribution of this signal is close to normal, in general. It can be shown that the plateau range in the histogram slope plot becomes more broad. It means that higher δ_max (a) and higher P_max (b) correspond to the same deflection error e in histogram slope plots. Thus, higher percent of distances is shorter than assumed δ_max which reduces calculation time.

There is the second gain of preliminary normalisation, as well. In the case of non-normalized signals the histogram convolutions must be calculated individually up to the highest embedding dimension used in calculations which takes a bit time. This is, however, unnecessary in the case of normalized signals because one can use the preliminary prepared histogram convolutions which do not change for every signals possessing normal distribution.

In the current section, I would like to show the interesting example of the signal possessing two plateaus in the $dw$=fn(W) relation. Let us consider the 6-dimensional system described by the equations Eq. 9 and presented in Figures 9 and 10:

The algorithm was applied to the variable y of this system. The results are presented in Figures 11 and 12. Two plateaus are visible in these Figures which make the interpretation difficult. The left plateau is visible at W ≈ 400. Let us observe that it corresponds to the elementary quasi-cycle of the attractor which is visible especially in variables z, $v$ and u (more than 40 quasi-cycles are observed in the period of 20000 samples). On the other hand, the long term plateau at W ≈ 1600 corresponds to about 4 elementary quasi-cycles in $w$⁠, $v$ and u or to a larger pack of shorter quasi-cycles in x, y and z. It can be interpreted that the second plateau determines the range of the higher regularity when compared to lower regularity observed at longer time scales.

This example shows that the problem of the estimation of the complexity is still open for new analyses and solutions, in general. The estimation of the correlation dimension is not a trivial task because it can be in many highly complex system time scale dependent. This example suggests the whole relation $dw$=fn(W) to be in some cases the descriptor of the system complexity rather than searching for the individual value of the correlation dimension of the signal.

The paper has presented some problems concerning the estimation of correlation dimension of HDS. The algorithm has been developed in order to calculate most precisely the correlation dimension using relative short time and to make the method acceptable for everyday use. The use and interpretation must be, however, careful for high dimensional signals, especially that possessing the correlation dimension higher than d > 6. The use of histogram convolution plots for the purpose of estimation of δ_max and $Pmaxo$ makes it possible to estimate the approximate error of calculations in the case of very high $dw$ > 10 which makes it possible to estimate it in the relative short time.

HDS

high dimensional signal,

d

correlation dimension of the signal,

$dw$

correlation dimension estimated using given W,

Embedding parameters

N

length of the analyzed signal,

L

lag or delay time,

m

embedding dimension,

W

window width, length of the extract of the signal creating one embedded point, W=(m − 1)·L,

n

number of embedded points,

A

matrix with embedded points [n×m],

Modified Takens-Ellner-Michalak algorithm parameters

LSA

“lower slope area” in the relation $dw$=fn(W), it corresponds to the plateau area in this function that is the feature of low-dimensional signals,

k =n/2

number of distances between pairs of points for 1 repetition of RJP,

δ

sorted vector with distances between embedded points,

I_min/I_max

index into vector δ pointing to the shortest/longest distance taken to d estimation,

D = I_max −I_min

number of distances used for d estimation,

D_a

assumed number of distances D to estimate d with assumed precision,

P_min/P_max

relative index into δ ranging between 0 and 1,

$Pmaxo$

approximate optimum P_max corresponding to the end of plateau in $dw$=fn(P_max) relation (based on assumed e),

δ_max

distance between embedded points corresponding to the end of plateau in $dw$=fn(P_max) relation (based on assumed e),

e

deflection error [%] in the histogram convolution plot. Approximate estimator of the error of d. $Pmaxo$ and δ_max can be approximately estimated for the assumed e.

W_opt

optimal W representing the point of minimum slope in the relation d=fn(W),

d_opt

correlation dimension calculated by providing W_opt to the 4-th (or other) degree polynomial fitted to $dw$=fn(W),

This section describes the Matlab function with algorithm implementation. The HDS-toolkit can be free downloaded from my private website: www.drmichalak.pl/chaos/eng/. The main function is “dhds.m” which can be just called as: d = dhds(signal) with default settings. However, a lot of parameters can be set individually depending on the measured correlation dimension, desired accuracy, desired computation time, memory use etc. Also, the intermediate parameters can be returned as output parameters, if necessary. Below, these parameters are described and its use is discussed.

The general call of this function is:[dcopt, Wopt, dcW, dcM, dcMc, WV, Pmax, other ] = dhds(series, quick_precise, MaxIterationTime, Da, DeflErr, norm, NW_Wmax_WV, Pmin, dcW, fitting, figures)

Use the empty string: ” (two apostrophes) to set the given parameter to be default and to define next parameters (e.g.):[dcopt, Wopt,dcW, dcM, dcMc, WV, Pmax, other] = dhds(series, ”, ”, ”, ”, 0) (normalization set to 0).

Let us analyze the input parameters:

• (1)

  quick_precise: this parameter makes it possible to set just one of 3 options of speed/precision: 1 or “quick” - quick, low precision, 2 or “medium” - medium (default), 3 or “precise” - slow, high precision. Use 1 for preliminary estimation and 3 for final estimation. Change other parameters if the results are not satisfactory.

• (2)

  MaxIterationTime: this parameter makes it possible to schedule the calculations with respect to its duration. It denotes the maximal duration (in seconds) of every W-iteration. There is more difficult to find distances that are shorter than assumed δ_max when the estimated correlation dimension for the given W increases and the calculations become significantly longer. The solution is to increase the accepted deflection error e which determines higher δ_max. The estimation of d becomes possible, however, it is performed with a bit lower precision due to higher e. After MaxIterationTime is reached without finding Da distances, the algorithm takes decision: (a) if more than 50% of D_a is reached until now, the calculation follows up to reaching D_a, however, e is increased for next iteration; if less than 50% distances is reached, W-iteration is repeated with higher e. e increases by 1% by default and cannot excess 30%. The default MaxIterationTime is 60 s, 200 s or 1000 s depending on quick_precise value. It is possible, also, to change the HEADER’s variable “UserAction.” If “UserAction”=0 than the calculations are repeated using increased DeflErr. If “UserAction”=1, the algorithm stops and asks for the user’s decision: “Follow” using current e or “Break and repeat using increased e.”

• (3)

  Da: assumed nbr of distances being shorter than the estimated distance δ_max. δ_max is estimated based on assumed deflection error e and used embedding dimension m. When omitted, default value of D_a depends on the used quick_precise parameter: D_a = 2000 or 5000 or 10000.

• (4)

  DeflErr(1:2): the vector with 2 elements. The first element describes the assumed deflection error e for the first iterations of calculations to estimate current δ_max. The default value is 3%. It is increased automatically when the correlation dimension for the given W increases and it is difficult to find assumed number of distances Da in the assumed MaxIterationTime. Increasing e reduces the precision of calculations. It must be pointed, however, that to small e may shift the results into “noise” region of d=fn(P_max) relation. You can use higher values initial values e.g. 5-10% if the expected correlation dimension d is very high and/or the signal is short and/or you want obtain robust results in very quick time. Use lower values if the correlation dimension of the signal is rather low, the signal is noise-free and you would like reduce the error connected with e-correction. High values accelerate strongly the calculations, especially if $dw$ is high, however, “e-correction error” occurs.

  The second element of this input parameter DeflErr(2) = e_cut represents the deflection error for which the δ_cut is estimated. It must be equal or greater than e. δ_cut determines the cut-off distance and thus the length of the actual vector δ in the memory space and the duration of its sorting, after D_a is reached. The distances shorter than δ_cut and longer than δ_max are the reserve if increase in e is necessary for the given W-iteration. The iteration has not to be repeated. It influences also the graphs being printed as a output of the function. If you use higher value, the relation d=fn(P_max) is broader. The acceptable values lie between 0 and 30. The unit is [%]. By default, it is 4%.

  Both these parameters are automatically increased if the assumed number of distances D_a shorter than δ_max is not reached in the assumed iteration time. The calculation for the given W-iteration is repeated using DeflErr increased by 1 making the estimation of $dw$ possible, however, with a bit lower precision.

• (5)

  norm: preliminary normalisation of the signal. Acceptable values are 0 and 1. If 1 (default) - the series is converted to possess a normal distribution (0, 1). Normalisation can reduce the calculation time, because the higher ratio of distances is shorter than assumed δ_max and P_max is higher for the assumed e. Be careful if the histogram of the signal is strongly non-gaussian, and especially, if the histogram possesses the maximum/a on its borders. Normalisation makes it also possible to use preliminary estimated histogram slopes (stored in deltamaxM_PmaxM_normal.mat file which is attached to the HDS-toolkit). They have not to be estimated individually for the given signal.

• (6)

  NW_Wmax_WV: Here, you can set the vector WV with your Window Width values used for calculations. You can just set the vector with the assumed W values. E.g. you want to use the same W values for different similar signals. If the number of elements is less than 4 the algorithm finds default values depending on the number of elements:

  • 1-element number: it is interpreted as NW

  • 2-element vector it is interpreted as: [NW, W_max]

  • 3-element vector: it is interpreted as: [NW, W_max, W_maxdensity]

    NW is the number of elements in WV vector (default is 40), W_max is the maximal value in WV (default is 10τ (autocorelation time of the signal)) and W_maxdensity denotes the value of W for which W values are chosen more dense. If you expect the range of plateau in $dw$=fn(W) then set W_maxdensity to be the center of the plateau in order to increase a bit the precision of polynomial fitting (default is 2.5τ).

  • if [0,W_max] - $dw$’s are calculated for increasing W’s calculated using a rule:

    W_next = W_previous * W_incr…up to W_max value starting from W=5samples. Use this option as the most universal and independent on τ. W_incr is equal to 1.1-1.3 depending on quick_precise value.

  • if “fntau” or “tau” or ”: (DEFAULT) sets 40 values W from 0.1τ to 10τ due to estimated autocorrelation time of the signal, unequally distributed with the maximum density for about 2.5τ.

    The distribution of default W-values can be modified in this method by changing the default parameters in the HEADER of the function: NW=40, defaultWmax_over_tau=10 and defaultWdense_over_tau=2.5.

• (7)

  P_min: This parameter determines the method of the choice of P_min parameter. This parameter does not influence strongly the result of final d estimation if P_max is high enough. As a rule, it is necessary to exclude the initial part of the vector δ with the very short distances which do not create a regular plateau. In the case of noise contaminated signal, P_min represents the range of the strong influence of noise.

  The interpretation of this parameter depends on its value: If 0 < P_min < 0.1 than it is treated as P_min. You can set this value if you have observed the begin of the influence of noise on the graph d=fn(P_max, W) in the preliminary run of the algorithm.If P_min > 0.1 - it determines the deflection error (in %) for which P_min should be estimated.If P_min=0 than no P_min-cut-off is used and d is estimated using all the distances shorter than δ_max (I_min=1). Use this value if you want see on the graph all the relation d=fn(P_max, W) with the noise contamination region.

  The default value for P_min is -1 which determines the automatic approximate determination of P_min, separately for every W iteration. P_min is chosen to be the smallest P_max for which the estimated d lies in the range 66-150% of the preliminary/roughly estimated d. Look for the variable PminCut=1.5 in the HEADER of the function to modify this range.

• (8)

  dcW: This parameter determines the method of the estimation $dw$ for the current W iteration from the relation d=fn(P_max).

  The question arises, which fragment of this relation should be the base for the estimation of $dw$? It is observed that the relation d=fn(P_max) possesses lacunarities in the plateau region. Thus, in order to increase the precision, a number of points from this relation should be taken to estimate $dw$ as a mean over these points. The parameter dcW determines the the choice of $Imaxdown−Imaxup$ range for these points.

  1. - range is maximally broad: $Imaxdown$ = $Iminopt$+200, $Imaxup$ = I_cut,

  2. - (default) range is broad: $Imaxdown$ = 0.05*$ImaxoptImaxup$ = 1.5*$Imaxopt$⁠,

  3. - range is medium broad: $Imaxdown$ = 0.12*$Imaxopt$⁠, $Imaxup$ = 1.2*$Imaxopt$⁠,

  4. - range is narrow: $Imaxdown$ = 0.3*$Imaxopt$⁠, $Imaxup$ = 1*$Imaxopt$⁠,

  5. - use 1 point only: I_max = $Imaxopt$⁠.

• (9)

  fitting: This parameter determines the method of polynomial fitting applied to the relation d=fn(W).if 0: look for optimal fitting over different polynomials and self-defined functions in the function “optimal_fitting.m” of the HDS toolkit. Here, you can define other functions than polynomials to find optimal fitting. The procedure works on symbolic expressions and solving differential equations, thus be careful with defining variables being optimized.if [4,5,6] (or higher): define the degrees of the polynomials to be fitted.if [4]: there is only 1 possible minimum in the 3-degree 1st derivative. You reach unambiguous result: 0 or 1 plateau and d_opt estimated.if [5] or [6]: there are 0,1 or 2 possible minima in the 4th- and 5th-degree 1st derivatives.if [7] or [8]: there are 0, 1, 2 or 3 possible minima in the 6th- and 7th-degree 1st derivatives.if [4,5] [4 6] [4 5 9] [4 6 8 10]: you can use every combination of the degrees to be used. The graphs and results will be printed for every polynomial used. Increase NW when using high degree polynomials to increase fitting precision. DEFAULT is [4 6 8].

• (10)

  figures:

  1. - turns off plotting figures

  2. - figures are plotted once after the last W-iteration.

  3. - turns on plotting intermediate figures after every W-iteration.

• I.

  dopt{i, j}(k): final estimators of correlation dimension

  • i=1 -

    results without e-correction, i=2 - results with e-correction

  • j -

    the number of the fitted polynomial, index into “fitting” vector

  • k -

    the number of plateaus for the given polynomial {i, j}. In case of higher degree polynomials, the number of the plateaus may be higher than 1.

• II.

  W_opt{i, j}(k): W’s corresponding to d_opt. i, j, k as above.

• III.

  dcW(1, :) - vector with $dw$’s estimated for subsequent W’s without e-correction. dcW(2,:) - with e-correction.

• IV.

  dcM - matrix (sizePmaxV, sizeWV) with d estimated without e-correction for the individual pairs P_max, W. Used for printing d=fn(P_max, W) relation.

• V.

  dcM_c - the same but with e-correction

• VI.

  WV - vector with W’s for which $dw$ were estimated to build the relation $dw$ = fn(W)

• VII.

  PmaxV - defined in the HEADER of the function as: PmaxV_default = 10^{(−8:0.01:−0.01)}. Necessary for printing graphs.

• VIII.

  other - structure keeping many other intermediate parameters. For details - see file dhds.m, lines 1081-1109.