On the Automatic Parameter Selection for Permutation Entropy

05/15/2019 ∙ by Audun Myers, et al. ∙ Michigan State University 0

Permutation Entropy (PE) has been shown to be a useful tool for time series analysis due to its low computational cost and noise robustness. This has drawn for its successful application in many fields. Some of these include damage detection, disease forecasting, and financial volatility analysis. However, to successfully use PE, an accurate selection of two parameters is needed: the permutation dimension n and embedding delay τ. These parameters are often suggested by experts based on a heuristic or by a trial and error approach. unfortunately, both of these methods can be time-consuming and lead to inaccurate results. To help combat this issue, in this paper we investigate multiple schemes for automatically selecting these parameters with only the corresponding time series as the input. Specifically, we develop a frequency-domain approach based on the least median of squares and the Fourier spectrum, as well as extend two existing methods: Permutation Auto-Mutual Information (PAMI) and Multi-scale Permutation Entropy (MPE) for determining τ. We then compare our methods as well as current methods in the literature for obtaining both τ and n against expert-suggested values in published works. We show that the success of any method in automatically generating the correct PE parameters depends on the category of the studied system. Specifically, for the delay parameter τ, we show that our frequency approach provides accurate suggestions for periodic systems, nonlinear difference equations, and ECG/EEG data, while the mutual information function computed using adaptive partitions provides the most accurate results for chaotic differential equations. For the permutation dimension n, both False Nearest Neighbors and MPE provide accurate values for n for most of the systems with n = 5 being suitable in most cases.



There are no comments yet.


page 1

page 2

page 3

page 4

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Permutation Entropy (PE) has it origins in information entropy, which is a tool to quantify the uncertainty in an information-based system. Information entropy was first introduced by Shannon [38]

in 1948 as Shannon Entropy. Specifically, Shannon entropy measures the uncertainty in future data given the probability distribution of the data types in the original, finite dataset. Shannon entropy is calculated as

, where represents a data type, and is the probability of that data type. In recent years information entropy has been heavily applied to the time series of dynamic systems. To better accommodate these applications, several new variations of information entropy have been proposed. Some of these include approximate entropy [32], sample entropy [35], and permutation entropy (PE) [2]. These methods measure the predictability of a sequence through the entropy of the relative data types. However, PE considers the ordinal position of the data (permutations), which have been shown to be very effective for analyzing the dynamic state of a time series [8, 29, 26]. PE is quantified in a similar fashion to Shannon entropy with only a change in the data type to permutations (see Fig. 2), which we symbolically represent as . PE has two parameters: the permutation dimension and embedding delay , which are used when selecting the permutation size and spacing, respectively. In practice, the results of PE have been shown to be highly dependent on these parameters [23, 39, 36]. Unfortunately, there is no single, accurate approach for selecting them. This introduces the motivation for this paper: investigate automatic methods for selecting both PE parameters. There are currently three main methods for selecting PE parameters: (1) use parameters suggested by experts for a specific application, (2) use trial and error to find suitable parameters, or (3) use methods developed for phase space reconstruction.

Many domain scientist who apply permutation entropy make general suggestions for and [44, 14]. However, these suggested parameters can result in inaccurate results. As an example of this shortcoming, Popov et al. [33] applied PE to EEG data and showed that the sampling frequency influences the selection of . As for the dimension , there are only general suggestions on how to choose its value. For example, it has been shown that the appropriate permutation dimension lies in the range for the vast majority of applications [36]. Additionally, Bandt and Pompe [2] suggest that , where is the length of the time series. However, these general outlines for the selection of (and ) do not allow for an application specific suggestion for .

If we assume that suitable PE parameters correspond to optimal phase space reconstruction parameters, then a common approach for selecting and

is to implement one of the existing methods for estimating the optimal Takens’ embedding 

[41] parameters. Based on this practice, some of the common methods for determining include the mutual information function approach [15], the first folding time of the autocorrelation function [16, 5], and phase space methods [7]. Additionally, some common phase space reconstruction methods for determining include box-counting [4], correlation exponent method [16], and false nearest neighbors [18]. However, although the parameters in PE have similar names to their delay reconstruction counterpart, there are innate differences between ordinal patterns and phase space reconstruction which can lead to inaccurate or values.

1.1 Our contribution

The problem we address in the this paper is the following: Given a sufficiently sampled/oversampled, noisy time series , how can we reliably and systematically define appropriate dimension and time delay values for computing the corresponding PE?

Our first contribution towards answering this question is detailed in Section 2 and it is related to the automatic selection of the time delay . Specifically, assuming that is contaminated by Gaussian measurement noise, in Section 2.1

we combine the Least Median of Squares (LMS) approach for outliers detection with Fourier transforms theorem to derive a formula for the maximum significant frequency in the Fourier spectrum. This formula allows obtaining a cutoff value where the only input, besides the signal, is a desired percentile from the Probability Density Function (PDF) of the Fourier spectrum. Once this value is obtained, Nyquist’s sampling theorem is used to compute an appropriate


The second contribution is through an approach that we develop in Section 2.2, which uses Multi-scale Permutation Entropy (MPE) for finding . We show how MPE can be used to find the main period of oscillation for a time series derived from a periodic system. Building upon this, we show how the method can be extended to find for a chaotic time series by using the first maxima in the MPE. Specifically, we suggest using the delay at the first maxima for PE calculations as it satisfies the Nyquist’s sampling theorem.

Our third contribution to the automatic selection of is through the analysis of Permutation Auto-Mutual Information [24] (PAMI). PAMI is an existing method for measuring the mutual information of permutations. However, in this paper we tailor this method to specifically select for PE. We found that PAMI reaches a first minima at a specific window size that is approximately independent of . We then use this window size to find a suitable based on the common selection of . We also conclude that finding the window size using is the least computationally expensive method when using PAMI.

Our final contribution towards answering the posited question is our evaluation of the ability of existing tools for computing an embedding dimension to provide an appropriate value for the PE parameter . We compare values computed from False Nearest Neighbors or (FNN—Section 3.1), Singular Spectrum Analysis (SSA—Section 3.2), and Multi-scale Permutation Entropy (MPE—Section 2.2). While we use existing methods for performing the FNN and the SSA analyses, for the MPE-based approach, we use a criteria established in prior works [36], which requires finding first, for selecting . However, for the latter, we further automate the selection process by using the procedure that we developed in Section 2.2 for finding first. The resulting and values from all the methods are then listed in Tables 2 and 3 and compared to values suggested in the literature for a variety of examples.

This paper is organized as follows. In Section 1.2 we provide a simple PE calculation example to illustrate the main concepts of PE and to introduce the influence of the two parameters and . We then go into detail on some existing methods for selecting both and . Specifically, in Section 2 we provide a detailed explanation for selecting using existing, automatic methods such as autocorrelation in Section 2.3 and Mutual Information (MI) in Section 2.4. Additionally, we modified and develop/tailor methods to automatically select . These methods include a frequency approach in Section 2.1, Multi-scale Permutation Entropy (MPE) in Section 2.2, and PAMI in Section 2.5. In Section 3 we expand on the process for selecting using False Nearest Neighbors (FNN) in Section 3.1 and Singular Spectrum Analysis in Section 3.2. In Section 3.3, we explain our algorithm for automatically selecting using MPE. After introducing each method, in Section 4 we contrast all of these methods and make conclusions on their viability by comparing the resulting parameters to those suggested by PE experts.

1.2 Permutation Entropy Example

Bandt and Pompe [2] defined PE according to


where is the probability of a permutation and is the permutation entropy for dimension , which has units of bits when the logarithm is of base 2. The permutation entropy parameters and are used when selecting the motif size, with as the time difference between two consecutive points in a uniformly sub-sampled time series and as the permutation length or motif dimension. To form a permutation, begin with with an element of the of the series . Using this element, the dimension , and delay

, define the vector

. The corresponding permutation of this vector is determined using its ordinal pattern. For example, consider the third degree permutation shown in Fig. 1. The permutation type, which categorizes the permutation, is found by first ordering the values of the permutation smallest to largest, and then accounting for the order received. For the given permutation in Fig. 1, the resulting permutation is categorized as the sequence , which is one of possible permutations for a dimension , see Fig. 2.

Figure 1: Sample permutation formation for and .
Figure 2: All possible permutation Configurations for n = 3.

To demonstrate how PE is calculated, consider another sequence with PE parameters and . The left side of Fig 3 shows how the sequence can be broken down into the following permutations: two , one , two , and one for a total of 6 permutations. This makes each permutation type have a probability out of 6. The permutation distribution can be visually understood by illustrating the probabilities of each permutation as separate bins. To accomplish this, the right side of Fig. 3 shows the abundance of each permutation.

Figure 3: Permutations 1 through 6 shown for example sequence (left) with and and the relative abundance of each permutation (right). Permutation 1 corresponds to a (0, 1, 2), permutation 2 is of type (0, 1, 2), permutation 3 is of type (1, 2, 0), permutation 4 is of type (1, 0, 2), permutation 5 is of type (1, 2, 0), and permutation 6 is of type (2, 1, 0).

Applying Eq. (1) to the probabilities of each permutation for our example sequence yields

We can normalize PE using the possible PE value, which occurs when all possible permutations are equiprobable according to . The resulting normalized PE is


Continuing with our example series , Eq. (2) gives a normalized PE .

2 Embedding Delay Parameter Selection

The delay embedding parameter is used to uniformly subsample the original time series. To elaborate, consider the time series . By applying the delay , a new subsampled series is defined as . In order to obtain a stable and automatic method for estimating an optimal value for we investigate: a novel frequency-based analysis that we describe in Section 2.1, Multi-scale Permutation Entropy (MPE) (Section 2.2), autocorrelation (Section 2.3), and Mutual Information function (MI) (Section 2.4). We recognize, but do not investigate, some other methods for finding such as diffusion maps [3] and phase space expansion [7].

2.1 Frequency Approach: Fourier Spectrum Noise Cutoff using -D Least Median of Squares

In this section, we develop a method for finding the noise floor in the Fourier spectrum using Least Median of Squares (LMS) [25]

. We then use the noise floor to find the maximum significant frequency of a signal contaminated with additive Gaussian white noise (GWN). Our method is based on finding the maximum significant frequency in the Fourier spectrum and the Nyquist sampling frequency criteria. To motivate the development of this approach, we begin by working with the frequency criteria developed by Melosik and Marszalek 

[27], which agrees with Nyquist sampling theorem [20], for choosing a suitable sampling frequency :


where is the maximum significant frequency in the signal. Melosik and Marszalek [27] showed that a sampling frequency within this range is appropriate for subsampling an oversampled signal, thus mitigating the effect of temporal correlations of neighboring points in densely sampled signals. They also showed that the time delay obtained using this criteria can also be used for delay reconstruction. However, the automatic identification of from an oversampled signal is not trivial. Melosik and Marszalek [27] selected a maximum significant frequency by inspecting the normalized Fourier spectrum and using a threshold cutoff of approximately for a noise-free chaotic Lorenz system. This made visually finding the maximum frequency significantly easier but did not provide guidance on how to algorithmically find . Further, Attempting to algorithmically adopt the approach in [27] resulted in large errors especially in the presence of a low signal to noise ratio. Therefore, this motivated the search for an automatic, data-driven approach for identifying the noise floor, which could then be used to find the maximum significant frequency. This led us to develop a method that is based on -D least median of squares applied to the Fourier spectrum, which we describe below. The assumptions inherent to our method are

  1. [nosep]

  2. The time series is not undersampled. The purpose of the methods is to determine a suitable delay parameter for subsampling the signal, which would be meaningless if the signal is undersampled.

  3. The time series has less than 50 outliers. By outliers, it is meant that the Fourier transform of the signal needs to have less than 50 of the points as significant peaks. This requirement stems from the limitations of the least median of squares regression.

  4. The noise in the signal is approximately Gaussian white noise; otherwise, the ensuing statistical analysis becomes inapplicable. Violating this assumption can yield false peak detections, which would lead to an incorrect delay parameter.

Using the noise floor determined from -D least median of squares, we find suitable cutoffs for obtaining of the signal and compute a suitable embedding delay according to


where we set , thus agreeing with the range in Eq. (3) and the Nyquist sampling criterion.

Figure 4 summarizes the frequency approach for

with the use of our 1-D LMS method for finding a noise floor in the Fourier spectrum. This process begins with computing the signal’s Fourier spectrum. We then fit an LMS regression line to noise in the Fourier spectrum to provide statistical information about the Probability Distribution Function (PDF) of the noise. Next, we use the PDF to determine the Cumulative Distribution Function (CDF), which can be used to determine a meaningful noise cutoff in the Fourier spectrum. However, if noise is present, it must be approximately GWN for this method to hold statistical significance. We then use this cutoff for separating the highest significant frequency in the Fourier spectrum

, which is used to find a suitable embedding delay based on the frequency criteria in Eq. (4). In the following paragraphs we review our use of the LMS and the derivation of the PDF of the Fourier spectrum of GWN. We then show how to combine the LMS method with the resulting PDF expression to find a suitable noise floor cutoff and the corresponding maximum significant frequency.

Figure 4: Overview of of our frequency domain approach for finding the maximum significant frequency using LMS for a signal contaminated with GWN.
Least Median of Squares:

LMS is a robust regression technique used when up to 50 of the data is corrupted by outliers. The algorithm for LMS regression was initially published in 1984 by Massart et al. [25]. In comparison to the widely used least sum of squares (LS) algorithm, the LMS replaces the sum for the median, which makes LMS resilient to outliers. The difference between LS and LMS can be described by


where is the residual. Figure 5 shows an example application of the linear LMS regression.

Figure 5:

LMS linear regression with 45

outliers. Results match those found in [25].

Specifically, this figure shows 110 data points drawn from the line

with added GWN of zero mean and 0.1 standard deviation. The data is corrupted with 90 outliers centered around

with a normal distribution of 1.0 along

and 0.6 along . Figure 5 shows that the linear regression results closely match the actual trend line with the fitted line being in comparison to the actual . These results agree well with the ones found in [25].

PDF and CDF of the magnitude of the Fast Fourier Transform of GWN:

This section reviews the probability distribution function (PDF) and cumulative density function (CDF) for the Fourier Transform (FT) of GWN. Additionally, this section derives the location of the theoretical maximum of the PDF. The distribution of the FT of GWN is described by [34]


where is the magnitude of the FT of GWN, is the probability density function of , is the standard deviation of the GWN, and is the window energy or number of discrete transforms taken during the FT. By setting the first derivative of with respect to equal to zero, the theoretical maximum of the PDF is


We calculate the CDF corresponding to the PDF described in Eq. (7) by combining the PDF in Eq. 6 with the CDF for a Rayleigh distribution in [30]. This results in


where is the cumulative probability of .

Finding the Noise Floor:

Our approach for finding the noise floor combines LMS with Eqs. (6) and (7). Specifically, we utilize LMS to obtain a 1-D fit of the Fast Fourier Transform (FFT) of the signal, which results in an approximate value of , which is at the maximum of . Using from the LMS fit, we then find the standard deviation of the distribution from Eq. 7, which is used to find a cutoff based on a set cumulative probability in Eq. (8).

We begin by showing the accuracy of the LMS fit for finding . Our example uses GWN with a mean of zero and standard deviation of 0.035 with 1000 data points. Taking the FFT of the GWN ( see Fig. 6A) results in the distribution shown in Fig. 6B. The distribution shows a 1-D LMS fit of 8.215 compared to the theoretical maximum of the PDF from Eq. 7 of 7.826, which is approximately 4.67 greater. This shows that the 1-D LMS fit accurately locates . Additionally, the theoretical shape of the PDF in Fig. 6B is shown to be very similar to the actual distribution.

Figure 6: (A) FFT of GWN with 0.035 standard deviation and zero mean with the location of the theoretical maximum of the PDF and one-dimensional LMS regression value. (B) Distribution of GWN in the Fourier Spectrum with overlapped theoretical PDF and location of the theoretical maximum of the PDF and one-dimensional LMS regression value.

Next, our approach utilizes Eq. (8) and derived from Eq. 7 for finding the cutoff value . The for a desired cumulative probability is found by solving Eq. (8) for as


In order to make robust to normalization and scaling of the FFT, we define the ratio between the suggested cutoff from Eq. (9) and the maximum of the PDF from Eq. (7) as

Example Cutoff:

An example of how Eqs. (7) and (9) are used is shown in Fig. 7, where the maximum of the PDF and the cutoff for are marked in Fig. 7a and 7b, respectively. For this example, we find the ratio C to be approximately 3.03 for a 99 probability. In addition, we suggested a cutoff ratio to be used for signals with less than data points. This yields an expected probability of approximately for a point in the FFT of the GWN attaining a magnitude greater than . Alternatively, Eq. (10) can be used to calculate a different value of C based on the desired probability and length of the signal.

Figure 7: (a) Theoretical PDF for GWN. (b) CDF for GWN with an example cutoff at the 99 .

2.2 Multi-scale Permutation Entropy for Embedding Delay

In this section, we develop a method based on Multi-scale Permutation Entropy (MPE) to find the periodicity of a signal, which is then used to find a suitable delay parameter. MPE is a method of applying permutation entropy over a range of delays, which was first introduced by Costa et al. [10] for analyzing physiological time series. Applying MPE results in a plot of permutation entropy versus delay, which reveals several useful trends. More specifically, Zunino et al. [45] showed how the first maxima in the MPE plot is developed when matches the characteristic time delay . Furthermore, the periodicity can be captured by the first dip in the MPE plot as shown in Fig. 8 at the location .

Figure 8: Resulting MPE plot for 2P periodic time series with example embedding vector sizes , , and .

Figure 8 shows embedding delays , , and calculated as as well as their corresponding locations on a normalized MPE plot. This MPE plot shows that the normalized MPE reaches its first maximum when the delay is roughly , which corresponds to to approximately an even distribution of motifs. A second observation, as mentioned previously, is that at (or the first dip in the MPE plot) there is a resonance or aliasing effect, which can be used to determine the period of the time series. This is based on the embedding delay size at causing the embedding vector size to be approximately half the of the periodicity P, which can be expressed as


where is the delay that causes aliasing, is the main period of oscillation, is the main frequency of the time series corresponding to , and is the sampling frequency. The reason for the dip in the permutation entropy (PE) when the condition from Eq. (11) is met is caused from an aliasing effect, which induces more regularity in the motif distribution.

Calculating Embedding Delay from MPE

We use the criteria of Melosik and Marszalek [27] to determine a suitable delay from the location of the first dip at . Their criteria states that the sampling frequency must fall within the range shown in Eq. (3). This range led to Eq. (4), which is used to calculate . However, for MPE, we substitute and in Eq. (3) with from Eq. (11) and . These substitutions allow Eq. (4) to reduce to


where . These simplifications show that is only dependent on the delay which causes resonance when applying MPE.

In the following subsections we will further explain the method, trends, and limitations for using MPE to calculate . Specifically, Section 2.2.1 explains in detail the trends in the MPE plot we use, Section 2.2.2 shows the MPE plot for an example chaotic time series from the Lorenz system, Section 2.2.3 investigates the robustness of the method to noise, and Section A details our algorithm (Algorithm 1) for finding using MPE.

2.2.1 MPE Regions

Riedl et al. [36] showed that the MPE plot can be separated into three distinct regions as described below and shown in Fig. 9.

Figure 9: The three regions of the MPE plot for a periodic signal: (A) redundant, (B) resonant, and (C) irrelevant.
Region A

Region A shows a gradual increase in the permutation entropy until reaching a maxima at the transition between regions A and B. Oversampling or a low value of causes the motif distriubtion corresponding to the permutation entropy to be heavily weighted on just increasing or decreasing motifs (motifs (0,1,2) and (2,1,0) for from Fig. 2). This effect was coined as the “Redundancy Effect" by De Micco et al. [13], which means sufficiently low values of result in redundant motifs. However, as increases, the motif distribution becomes more equiprobable. Additionally, when the motif probability reaches a maximum equiprobability, the permutation entropy is at a maxima, which is the point of transitions from region A to B.

Region B

Region B shows a slight dip to the first minima. This reduction in permutation entropy is caused by the aliasing or resonance from the value of approaching half the main period length. At the transition from B to C, the resonance is reached, which provides information on the main frequency and period of the time series.

Region C

Region C has possible additional minima and maxima from additional alignment of the embedding vector with multiples of the main period. This region was referred to as the “Irrelevant Region" by De Micco et al. [13] due to effectively large values of forcing the delayed sampling frequency to fall below the Nyquist sampling rate as described by the lower bound in Eq. (3).

2.2.2 Example with Chaotic Time Series

Figure 10: MPE plot for the coordinate of the Lorenz system. Additionally, points in the MPE plot with their corresponding subsampled time series are shown for the redundant, resonant, and irrelevant regions as described in Section 2.2.1.

In Sections 2.2 and 2.2.1, we used a periodic time series to show and explain the regions developed in an MPE plot as well as an MPE-based method for determining a suitable embedding delay . In this section we further show the applicability of this approach to chaotic signals using the -coordinate of the Lorenz System as an example. We simulate the Lorenz equations


with a sampling rate of Hz and using the parameters , , and . This system was solved for seconds and only the last seconds from the time series are used. Figure 10 shows the result of applying MPE to the simulated Lorenz system.

Figure 10 shows similarities to Fig. 9 with a clear maxima at the boundary between regions A and B, albeit with no obvious minima. Therefore, a new distinct feature needs to be used to determine . We suggest using the first maxima to find because this delay is likely to fall within the region described by Eq. (12).

2.2.3 Effects of Noise

Figure 11: Region N is affected by noise in the MPE plot, and region S is unaffected.

We found that the main advantage of using MPE for determining the embedding delay is its robustness to noise. Noise on an MPE plot has minimal effects on regions B and C from Fig. 9, while only significantly affecting region A as shown in Fig 11. Furthermore, depending on the signal to noise ratio, there will only be an effect at the beginning of region A. Figure 11 shows the first region N where noise is affecting the permutation entropy. The effect of noise causes the MPE plot to start at a maxima and decrease to a local minima. When the time delay becomes large enough, the permutations are no longer influenced by the noise causing this minima. We found that the location of the minima is based on the condition


where is the average of the absolute value of the slope and is approximately the maximum amplitude of the noise, is the value of great enough to surpass the noise amplitude. We derived this condition from the need for, on average, . This shows that MPE is robust to noise as long as the noise amplitude does not exceed the amplitude of the signal.

2.3 Autocorrelation

Autocorrelation uses the correlation coefficient between the time series and its -lagged version to find a suitable delay for phase space reconstruction. This method was first introduced by Box et al. [5] in 1970. Typically, the autocorrelation function is computed as a function of and, as a rule of thumb, a suitable delay is found when the correlation between and reaches the first folding time, i.e., when  [17]. We will overview the two prominent correlation techniques that can be used when implementing an autocorrelation-based approach for finding . These methods include Pearson Correlation in section 2.3.1 and Spearman’s Correlation in section 2.3.2. After these two methods are introduced, we provide an example demonstrating the use of autocorrelation in section 2.3.3.

2.3.1 Pearson Correlation

The Pearson correlation coefficient measures the linear correlation of two time series and . Using these two data sets the correlation coefficient is calculated as


The possible values of represent the relationship between the two data sets, where represents a perfect positive linear correlation, represents no linear correlation, while represents a perfect negative linear correlation. However, Pearson correlation is limited because it only detects linear correlations. This limitation is somewhat alleviated by using Spearman’s Correlation which operates on the ordinal ranking of the two time series instead of their numeric values.

2.3.2 Spearman’s Correlation

Spearman’s correlation is also calculated using Eq. (15) with the substitution of and for their ordinal ranking. This substitution allows for detecting nonlinear correlation trends to be represented as long as the correlation is monotonic. To demonstrate the difference, Fig. 12 shows two sequences and calculated from with . Using this example, the Pearson correlation is calculated as , while Spearman’s ranked correlation yields . This result demonstrates how Spearman’s correlation coefficient accurately detects the non-linear, monotonic correlation between and whereas Pearson correlation may miss it.

Figure 12: A comparison between (left) unranked values and (right) ranked values for calculating correlation coefficients. Using the ranked and , Spearman’s correlation coefficient can be used to accurately reveal existing nonlinear monotonic correlations.

2.3.3 Autocorrelation Example

We can use the concept of correlation to select a delay by calculating the correlation coefficient using Eq. (15) between a time series and its -lagged version. As an example, take the time series , with having a sampling frequency of Hz. This results in a suggested delay at the first folding time using both Spearman’s and Pearson correlation. In section 4 we will implement Spearman’s version of autocorrelation to account for the possibility of non-linear correlations.

2.4 Mutual Information

Mutual information (MI) can be used to select the embedding delay based on a minimum in the joint probability between two sequences. The mutual information between two discrete sequences was first realized by Shannon et al. [37] as


where and are the two sequences, and are the probability of the element and separately, and is the joint probability of and . Fraser and Swinney [15] showed that for a chaotic time series the MI between the original sequence and and delayed version will decrease as increases until reaching a first minimum. At this minima, the delay allows for the individual data points to share a minimum amount of information, which indicates sufficiently separated data points. While this delay value was specifically developed for phase space reconstruction, it is also commonly used for the selection of the PE parameter . We would like to point out that, in general, there is no guarantee that local minima exist in the mutual information, which is a serious limitation for computing using this method. All MI methods can be applied to either ranked or unranked data, and we will investigate the resulting delay values using both of these scenarios. As an overview, we plan to investigate four methods for estimating for PE using MI. These methods include MI with equal-sized partitions in Section 2.4.1, adaptive partitions in Section 2.4.2, and two permutation-based MI estimation methods in Section 2.4.3. After presenting each method, In Section 2.4.4 the optimal will be chosen by comparing the values from each MI method to those suggested by PE experts for various systems.

2.4.1 MI using Equal-sized Partitions

For the calculation of MI, the joint and independent probabilities of the original and time lagged time series are needed. However, since is a discrete time series, we approximate these probabilities using bins, which segment the range of the series into discrete groups. The simplest method for approximating the probabilities using this discretization method is to use equal sized bins. However, the size of these bins is dependent on the number of bins . We investigated various methods for estimating an appropriate number of bins using the length of the time series . These methods include the common square-root choice , Sturge’s formula [40] , and Rice Rule [21] . After comparing each method using a variety of examples, we found that the use of Sturge’s formula provided the best results for selecting for PE using MI.

2.4.2 MI using Adaptive Partitions

Darbellay and Vajda [12] introduced a multistep, adaptive partitioning scheme to select appropriate binning sizes in the observation space formed by the plane and . Their method is often considered state-of-the-art for estimating the mutual information function [19]. In this approach, the bins are recursively created where in the first function call, the space of the signal and its -lagged version is divided into an equal number of

D bins. Then a A chi-squared test is used to test the null hypothesis that the data within the newly created bins are independent. Any segment that fails the test is further divided until the resulting sub-segments contain independent data (or a certain number of divisions is satisfied). Using this partitioning method, the MI is calculated using Eq. (


2.4.3 Kraskov MI

Kraskov et al. [19] developed a method for approximating the MI using entropy estimates using partition sizes based on -nearest neighbors. Specifically, the method begins by first calculating the MI using entropy [11] as


where is the Shannon entropy. Next, an approximation of with digamma functions is done, but the probability density of and still needs to be estimated. To do this, adaptive partitions using the -nearest neighbor are formed. Specifically Kraskov et al. develop two different partitioning methods with similar results. The first method uses the maximum Chebyshev distance to the nearest neighbor to form square bins as shown in Fig. 13-a, and the second method in Fig. 13-b uses rectangular partitions using the horizontal and vertical distances to the nearest neighbor .

Figure 13: Example showing two different partition methods for Mutual Information estimation using an using nearest neighbors adaptive partitioning.

To continue with the example shown in Fig. 13, the density probability is estimated using the strips formed from these bins. To highlight the difference, Fig. 13-a shows a horizontal strip of width encapsulating points (strip does not include the point ), while in Fig. 13-b only point is enclosed. Using these probability density approximations and the digamma function , MI between and can be estimated. Using the partitioning method shown in Fig. 13-a the MI is estimated as


Using the partitioning method shown in Fig. 13-b the MI is estimated as


For the results shown in Section 4, we use the nearest neighbor to generate the partitions.

Mutual Information
Kraskov et al.
Method 1
Kraskov et al.
Method 2
White Noise 1 3 3 1 1 [36]
Lorenz 13 9 9 9 10 [36]
Rossler 14 13 11 9 9 [42]
16 14 14 15 15 [36]
Mackey-Glass 7 8 7 7 1 to 700 [36]
Sine Wave 4 17 13 1 15 [42]
Logistic Map 5 8 11 5 1 to 5 [36]
Henon Map 12 15 13 8 1 to 5 [36]
ECG 22 16 9 8 1 to 4 [36]
EEG 6 5 5 5 1 to 3 [36]
Table 1: A comparison between the calculated and suggested values for the delay parameter for multiple MI approximation methods. The shaded cells highlight the methods that yielded the closest match to the suggested delay. The equal-sized partition method is described in Section 2.4.1, Kraskov et al. methods and in Section 2.4.3, and the adaptive partitioning approach in Section 2.4.2.

2.4.4 MI Methods Comparison

To determine the optimal MI approximation method for selecting for PE, Table 1 shows a comparison between the values computed from each of the MI methods in Sections and the corresponding values suggested by experts. The table shows that the adaptive partitioning method of Section 2.4.2 results in an accurate selection of for the majority of systems. Therefor, we will use the adaptive partitioning estimation method when making comparisons to other methods in Section 4.

2.5 Permutation Auto-mutual Information

As shown in Section 2.4, Mutual information (MI) is a useful method for selecting for phase space reconstruction. However, it does not account for permutations when selecting , which can lead to inaccuracies in computing the PE. To circumvent this issue, we develop a new method for selecting using Permutation Auto-Mutual Information (PAMI). PAMI was first introduced by Liang et al. [24] as a tool for detecting dynamic changes in brain activity. However, in this paper we are tailoring its application for the first time to the selection of the permutation entropy parameter . PAMI is defined as


where is the permutation entropy described in Eq. (1). We suggest an optimal delay for a given dimension when PAMI is at a minimum. This delay corresponds to minimum shared information between the time series and its time lagged version. By applying this method for the simple sinusoidal function described in Section A.2.3 we can form Fig. 14 with and .

Figure 14: PAMI results for the sinusoidal function described in Section A.2.3 with and . The figure shows an optimal window size .

As shown, the window size is approximately independent of the dimension , with an optimal window for the example. Through our analysis of the minimum PAMI as a function of the window size, we have developed a new method for selecting the optimal embedding window. However, we need the embedding dimension to suggest an optimal delay. To do this, we implement the common choice for ranging from for PE. To reduce the computational demand, we suggest using permutation dimensions to find an optimal window size. In addition to the reduced computational demand of using , we found that at the first minima. This also helps making this first minima even more simple.

3 Motif Dimension

The second parameter for permutation entropy that needs to be automatically identified is the embedding dimension . The methods for determining fall into one of two categories: (1) independently determining and , and (2) simultaneously determining and based on the width of the embedding window. For the first category, we investigate using the method of False Nearest Neighbors (FNN) [18] in Section 3.1, and Singular Spectrum Analysis (SSA)[6] in Section 3.2. For the second category, we contribute to the selection of by developing an automatic method using MPE from Section 3.3. This method combines the results for finding through MPE in Section2.2 with the work of Riedl et al. [36]. We recognize that our work does not include other commonly used methods for independently calculating such as box-counting [9], largest Lyapunov exponent [43], and Kolmogorov–Sinai entropy [31]. However, our comparisons include many of the most common ones.

3.1 False Nearest Neighbors

False Nearest Neighbors (FNN) is one of the most commonly used methods for determining the minimum embedding dimension for state space reconstruction. It was originally developed by Kennel et al. [18] as a geometric approach for determining the minimum needed embedding dimension. In this method the time series is repeatedly embedded into an a sequence of -dimensional Euclidean spaces for a range of increasing values of . The idea is that when the minimum embedding dimension is reached or , the distance between neighboring points does not significantly change as we keep increasing . In other words, the Euclidean distance between the point and its nearest neighbor minimally changes when the embedding dimension increases to . If the dimension is not sufficiently high, then the points are false neighbors if their pairwise distance significantly increases when incrementing . This ratio of change in the distance between nearest neighbors embedded in and is quantified using the ratio of false nearest neighbors


is compared to the tolerance threshold to distinguish false neighbors when . In this paper, we select as used by Kennel et al. [18]. By applying this threshold over all points, we can find the number of false neighbors as a percent FNN . If there is no noise in the system, should reach zero when a sufficient dimension is reached. However, with additive noise present, may never reach zero. Thus, it is commonly suggested to use a percent FNN cutoff for finding a sufficient dimension . We use the typically chosen cutoff , which is suitable for most applications when moderate noise is present.

3.2 Singular Spectrum Analysis

The singular spectrum analysis method was first introduced by Broomhead and King [6] as a tool to find the trends and most prominent periods of a time series. Leles et al. [22]

summarized the SSA procedure as (1) immersion, (2) Singular Value Decomposition (SVD), (3) grouping, and (4) diagonal averaging. Specifically, immersion embeds the time series into a dimension

to form a Hankel matrix, SVD factors all the matrices, grouping combines the matrices that are similar in structure, and diagonal averaging reconstructs the time-series using the combined matrices. The needed embedding dimension is determined from the SVD by calculating the ratio


of the sum of the th diagonal entries to the sum of the total diagonal entries . When exceeds , we consider the dimension to be high enough and set , which can then be used as the embedding dimension for permutation entropy.

3.3 Multi-scale Permutation Entropy

Riedl et al. [36] showed how MPE can be used to determine an embedding dimension . This method requires the embedding delay to be set to the length of the main period of the signal as shown in Section 2.2. The theory behind the method is based on normalizing the MPE according to


where is the PE normalized using the embedding dimension, and is the PE calculated from Eq. (1). Riedl et al. [36] determine the embedding dimension by incrementing to find the largest corresponding normalized PE with an embedding delay heuristically determined from the main period length. They concluded that the with the highest entropy accurately accounts for the needed complexity of the time series, and therefore suggests a suitable embedding dimension. Rield et al. [36] show how this method provides an accurate embedding dimension for the Van-der-Pol-oscillator, Lorenz system, and the logistic map. However, the method is not automatic due to the reliance on a heuristically chosen .

To make the process automatic, we introduce algorithm 1 based on Section 2.2 to automatically select the correct , which we then use in conjunction with Eq. (23) to find corresponding to the maximum . Additionally, we suggest scaling from to as we have not yet found a system requiring using this method.

4 Results and Discussion

To make conclusions about the described methods for determining and , we made comparisons to values suggested by experts. The majority of the suggested parameters are taken from the work of Riedl et al. [36]. Tables 2 and 3 show the calculated and suggested values for and , respectively.

Embedding Delay

Table 2 highlights the automatically computed which match relatively well with the expert-identified value for a variety of systems. These systems fall within several categories, which include the following: noise, chaotic differential equations, periodic systems, nonlinear difference equations, and medical data. The methods presented in Table 2 include PAMI from Section 2.5, MI calculated using adaptive partitioning from Section 2.4.2, Spearman’s Autocorrelation from Section 2.3, MPE from Section 2.2, and the frequency approach from Section 2.1. For the noise category we only investigated GWN, and all the methods accurately suggest an embedding delay. For the second category of chaotic differential equations, Mutual Information approximated using adaptive partitions accurately provided suitable delay values. For the third category, periodic systems, we only investigated a simple sinusoidal function. This resulted in both MPE and the Frequency approach providing accurate suggestions. For difference equations we found that PAMI, autocorrelation, MPE, and the frequency approach provide accurate suggestions for the delay. Finally, when testing each method on medical data with intrinsic noise, we found that the noise-robust frequency approach yielded the optimal parameter selection for . As a generalization of the results found, we suggest the use of MI with adaptive partitioning when selecting for chaotic differential equations. However, for periodic systems, nonlinear difference equations, or ECG/EEG data we suggest the use of the frequency approach that we developed in this paper.

Embedding Dimension

Table 3 shows parameters falling within the region specified by experts as highlighted in green. It can be seen that both MPE and FNN commonly had parameters within the range specified for all categories. However, SSA failed to provide a consistently suitable embedding dimension . This leads to the conclusion that either MPE or FNN are sufficient methods for determining the embedding dimension for the majority of the considered applications. Although, these results also show that the dimension works well for almost all applications.

Catagory System
MI using
Noise White Noise 1 1 1 1 1 1 [36]
Lorenz 5 to 9 9 15 17 6 10 [36]
Rossler 6 to 10 9 12 19 7 9 [42]
6 to 10 15 12 20 7 15 [36]
Mackey-Glass 2 to 4 7 5 8 3 1 to 700 [36]
Sine Wave 5 to 8 1 10 16 21 15 [42]
Logistic Map 1 5 1 1 1 1 to 5 [36]
Difference Eq.
Henon Map 1 8 1 1 1 1 to 5 [36]
ECG 1 to 2 8 21 13 2 1 to 4 [36]
EEG 2 to 4 5 4 4 1 1 to 3 [36]
Table 2: A comparison between the calculated and suggested values for the delay parameter . The shaded cells highlight the methods that yielded the closest match to the suggested delay. The following conditions or abbreviations were used in the table: the range under PAMI results is from using the range (), AP under MI is an abbreviation for adaptive partitioning, and AC is an abbreviation for autocorrelation.
Category System Method
5 4 23 3 to 6  [36]
Lorenz 5 3 4 5 to 7  [36]
Rössler 4 4 4 6  [42]
Rössler system
4 4 4 6 to 7  [36]
4 4 6 4 to 8  [36]
Sine Wave 3 4 2 4  [42]
Logistic Map 5 4 3 2 to 16  [36]
Henon Map 5 4 2 3 to 10  [36]
ECG 5 7 8 3 to 7  [36]
EEG 6 5 11 3 to 7  [36]
Table 3: A comparison between the calculated and suggested values for the embedding dimension . The shaded cells highlight the methods that yielded the closest match to the suggested dimension.

5 Conclusion

In this paper we demonstrated various methods for automatically determining the PE parameters and when supplied with a sufficiently sampled/oversampled time series. The goal is to find, in an automatic way, the most accurate method in comparison to expert suggested parameters. The methods we investigated for calculating include autocorrelation, mutual information, permutation-auto-mutual information, frequency analysis, and multi-scale permutation entropy. Additionally, the methods we investigated for determining the embedding dimension include false nearest neighbors, singular spectrum analysis, and multiscale permutation entropy.

Our first contribution was developing a new frequency approach analysis and extending two existing methods, PAMI, and MPE, to automatically determine . For the frequency approach, we developed an automatic algorithm for finding the maximum significant frequency using a cutoff greater than the noise floor. The noise floor was found using one dimensional least median of squares applied to the Fourier spectrum in conjunction with the theoretical probability distribution function for the Fourier transform of Gaussian white noise. For PAMI, we showed how the minimum in the PAMI can be used to find an optimal embedding window, which is approximately independent of the dimension . We then suggested using the range to find from the first minima in PAMI for . For MPE, we showed how it can be used to find the main period of oscillation from a periodic time series, which we then use to find . Additionally, we expanded upon this method by showing how the main period of oscillation can also be found for non-periodic time series, which we implemented into an automatic algorithm.

Our second contribution was implementing the automatic selection of using MPE to also find using MPE. We also collected and compared some of the most popular methods for obtaining including false nearest neighbors, and singular spectrum analysis. We applied these methods to various categories including difference equations, chaotic differential equations, periodic systems, EEG/ECG data, and Gaussian noise. We then compared the generated parameters to values suggested by experts to determine which methods consistently found accurate values for and . We found that SSA did not provide suitable values for . However, both FNN and MPE provided accurate values for for most of the systems. We also concluded that, for the majority applications, a permutation dimension is suitable. For determining , we showed that our frequency approach provided accurate suggestions for for periodic systems, nonlinear difference equations, and medical data, while the mutual information function computed using adaptive partitions provided the most accurate results for chaotic differential equations.


FAK acknowledges the support of the National Science Foundation under grants CMMI-1759823 and DMS-1759824.


  • [1] Ralph G Andrzejak, Klaus Lehnertz, Florian Mormann, Christoph Rieke, Peter David, and Christian E Elger. Indications of nonlinear deterministic and finite-dimensional structures in time series of brain electrical activity: Dependence on recording region and brain state. Physical Review E, 64(6):061907, 2001.
  • [2] Christoph Bandt and Bernd Pompe. Permutation entropy: a natural complexity measure for time series. Physical review letters, 88(17):174102, 2002.
  • [3] T. Berry, J. R. Cressman, Z. Gregurić-Ferenček, and T. Sauer. Time-scale separation from diffusion-mapped delay coordinates. SIAM Journal on Applied Dynamical Systems, 12(2):618–649, jan 2013.
  • [4] A Block, W Von Bloh, and HJ Schellnhuber. Efficient box-counting determination of generalized fractal dimensions. Physical Review A, 42(4):1869, 1990.
  • [5] George EP Box, Gwilym M Jenkins, Gregory C Reinsel, and Greta M Ljung. Time series analysis: forecasting and control. John Wiley & Sons, 2015.
  • [6] David S Broomhead and Gregory P King. Extracting qualitative dynamics from experimental data. Physica D: Nonlinear Phenomena, 20(2-3):217–236, 1986.
  • [7] Th Buzug and G Pfister. Optimal delay time and embedding dimension for delay-time coordinates by analysis of the global static and local dynamical behavior of strange attractors. Physical review A, 45(10):7073, 1992.
  • [8] Yinhe Cao, Wen-wen Tung, JB Gao, Vladimir A Protopopescu, and Lee M Hively. Detecting dynamical changes in time series using the permutation entropy. Physical review E, 70(4):046217, 2004.
  • [9] Septima Poinsette Clark. Estimating the fractal dimension of chaotic time series. Lincoln Laboratory Journal, 3(1), 1990.
  • [10] Madalena Costa, Ary L Goldberger, and C-K Peng. Multiscale entropy analysis of complex physiologic time series. Physical review letters, 89(6):068102, 2002.
  • [11] Thomas M Cover and Joy A Thomas. Elements of information theory. John Wiley & Sons, 2012.
  • [12] G.A. Darbellay and I. Vajda. Estimation of the information by an adaptive partitioning of the observation space. IEEE Transactions on Information Theory, 45(4):1315–1321, may 1999.
  • [13] Luciana De Micco, Juana Graciela Fernández, Hilda A Larrondo, Angelo Plastino, and Osvaldo A Rosso. Sampling period, statistical complexity, and chaotic attractors. Physica A: Statistical Mechanics and its Applications, 391(8):2564–2575, 2012.
  • [14] Birgit Frank, Bernd Pompe, Uwe Schneider, and Dirk Hoyer. Permutation entropy improves fetal behavioural state classification based on heart rate analysis from biomagnetic recordings in near term fetuses. Medical and Biological Engineering and Computing, 44(3):179, 2006.
  • [15] Andrew M Fraser and Harry L Swinney. Independent coordinates for strange attractors from mutual information. Physical review A, 33(2):1134, 1986.
  • [16] Peter Grassberger and Itamar Procaccia. Measuring the strangeness of strange attractors. Physica D: Nonlinear Phenomena, 9(1-2):189–208, 1983.
  • [17] Holger Kantz and Thomas Schreiber. Nonlinear time series analysis, volume 7. Cambridge university press, 2004.
  • [18] Matthew B Kennel, Reggie Brown, and Henry DI Abarbanel. Determining embedding dimension for phase-space reconstruction using a geometrical construction. Physical review A, 45(6):3403, 1992.
  • [19] Alexander Kraskov, Harald Stögbauer, and Peter Grassberger. Estimating mutual information. Physical Review E, 69(6), jun 2004.
  • [20] HJ Landau. Sampling, data transmission, and the nyquist rate. Proceedings of the IEEE, 55(10):1701–1706, 1967.
  • [21] David Lane, Joan Lu, Camille Peres, Emily Zitek, et al. Online statistics: An interactive multimedia course of study. Retrieved January, 29:2009, 2008.
  • [22] Michel CR Leles, João Pedro H Sansão, Leonardo A Mozelli, and Homero N Guimarães. Improving reconstruction of time-series based in singular spectrum analysis: A segmentation approach. Digital Signal Processing, 77:63–76, 2018.
  • [23] Duan Li, Zhenhu Liang, Yinghua Wang, Satoshi Hagihira, Jamie W. Sleigh, and Xiaoli Li. Parameter selection in permutation entropy for an electroencephalographic measure of isoflurane anesthetic drug effect. Journal of Clinical Monitoring and Computing, 27(2):113–123, dec 2012.
  • [24] Zhenhu Liang, Yinghua Wang, Gaoxiang Ouyang, Logan J Voss, Jamie W Sleigh, and Xiaoli Li. Permutation auto-mutual information of electroencephalogram in anesthesia. Journal of Neural Engineering, 10(2):026004, feb 2013.
  • [25] Desire L Massart, Leonard Kaufman, Peter J Rousseeuw, and Annick Leroy. Least median of squares: a robust method for outlier and model error detection in regression and calibration. Analytica Chimica Acta, 187:171–179, 1986.
  • [26] Michael McCullough, Michael Small, Thomas Stemler, and Herbert Ho-Ching Iu. Time lagged ordinal partition networks for capturing dynamics of continuous dynamical systems. Chaos: An Interdisciplinary Journal of Nonlinear Science, 25(5):053101, 2015.
  • [27] Michał Melosik and W Marszalek. On the 0/1 test for chaos in continuous systems. Bulletin of the Polish Academy of Sciences Technical Sciences, 64(3):521–528, 2016.
  • [28] George B Moody and Roger G Mark. The impact of the mit-bih arrhythmia database. IEEE Engineering in Medicine and Biology Magazine, 20(3):45–50, 2001.
  • [29] Audun Myers, Elizabeth Munch, and Firas A Khasawneh. Persistent homology of complex networks for dynamic state detection. arXiv preprint arXiv:1904.07403, 2019.
  • [30] Athanasios Papoulis and S Unnikrishna Pillai.

    Probability, random variables, and stochastic processes

    Tata McGraw-Hill Education, 2002.
  • [31] Yakov Borisovich Pesin. Characteristic lyapunov exponents and smooth ergodic theory. Uspekhi Matematicheskikh Nauk, 32(4):55–112, 1977.
  • [32] Steven M Pincus. Approximate entropy as a measure of system complexity. Proceedings of the National Academy of Sciences, 88(6):2297–2301, 1991.
  • [33] Anton Popov, Oleksii Avilov, and Oleksii Kanaykin. Permutation entropy of eeg signals for different sampling rate and time lag combinations. In Signal Processing Symposium (SPS), 2013, pages 1–4. IEEE, 2013.
  • [34] Mark A Richards. The discrete-time fourier transform and discrete fourier transform of windowed stationary white noise. Georgia Institute of Technology, Tech. Rep, 2013.
  • [35] Joshua S Richman and J Randall Moorman. Physiological time-series analysis using approximate entropy and sample entropy. American Journal of Physiology-Heart and Circulatory Physiology, 278(6):H2039–H2049, 2000.
  • [36] Müller Riedl, A Müller, and N Wessel. Practical considerations of permutation entropy. The European Physical Journal Special Topics, 222(2):249–262, 2013.
  • [37] Claude E Shannon, Warren Weaver, and Arthur W Burks. The mathematical theory of communication. 1951.
  • [38] Claude Elwood Shannon. A mathematical theory of communication. ACM SIGMOBILE mobile computing and communications review, 5(1):3–55, 2001.
  • [39] MATTHÄUS STANIEK and KLAUS LEHNERTZ. PARAMETER SELECTION FOR PERMUTATION ENTROPY MEASUREMENTS. International Journal of Bifurcation and Chaos, 17(10):3729–3733, oct 2007.
  • [40] Herbert A Sturges. The choice of a class interval. Journal of the american statistical association, 21(153):65–66, 1926.
  • [41] Floris Takens. Detecting strange attractors in turbulence. In Dynamical systems and turbulence, Warwick 1980, pages 366–381. Springer, 1981.
  • [42] Mei Tao, Kristina Poskuviene, Nizar Alkayem, Maosen Cao, and Minvydas Ragulskis. Permutation entropy based on non-uniform embedding. Entropy, 20(8):612, 2018.
  • [43] Alan Wolf, Jack B Swift, Harry L Swinney, and John A Vastano. Determining lyapunov exponents from a time series. Physica D: Nonlinear Phenomena, 16(3):285–317, 1985.
  • [44] Hong Zhang and Xuncheng Liu. Analysis of parameter selection for permutation entropy in logistic chaotic series. In Intelligent Transportation, Big Data & Smart City (ICITBS), 2018 International Conference on, pages 398–402. IEEE, 2018.
  • [45] Luciano Zunino, Miguel C Soriano, Ingo Fischer, Osvaldo A Rosso, and Claudio R Mirasso. Permutation-information-theory approach to unveil delay dynamics from time-series analysis. Physical Review E, 82(4):046212, 2010.

Appendix A Appendix

a.1 Dynamic System Models

a.1.1 Lorenz System

The Lorenz system used is defined as


The Lorenz system had a sampling rate of 100 Hz with parameters , , and . This system was solved for 100 seconds and the last 24 seconds were used.

a.2 Rössler System

The Rössler system used was defined as


with parameters of , , , which was solved over 400 seconds with a sampling rate of 10 Hz. Only the last 1500 data points of the x-solution were used in the analysis.

a.2.1 Bi-Directional Coupled Rössler System

The Bi-directional Rössler system is defined as


with , , and . This was solved for 4000 seconds with a sampling rate of 10 Hz. Only the last 400 seconds of the x-solution were used in the analysis.

a.2.2 Mackey-Glass Delayed Differential Equation

The Mackey-Glass Delayed Differential Equation is defined as


with , , , and . This was solved for 400 seconds with a sampling rate of 100 Hz. The solution was then downsampled to 5 Hz and only the last 1500 terms of the x-solution were used in the analysis.

a.2.3 Periodic Sinusoidal Function

The sinusoidal function is defined as


This was solved for 10 seconds with a sampling rate of 50 Hz.

a.2.4 EEG Data

The EEG signal was taken from andrzejak et al. [1]. Specifically, the first 2000 data points from the EEG data of a healthy patient from set A, file Z-093 was used.

a.2.5 ECG Data

The Electrocardoagram (ECG) data was taken from SciPy’s misc.electrocardiogram data set. This ECG data was originally provided by the MIT-BIH Arrhythmia Database [28]. We used data points 3000 to 4500 during normal sinus rhythm.

a.2.6 Logistic Map

The logistic map was generated as


with and . Equation 29 was solved for the first 500 data points.

a.2.7 Hénon Map

The Hénon map was solved as


where , , , and . This system was solved for the first 500 data points of the x-solution.

a.3 Algorithm

a.3.1 MPE delay Algorithm

In Section 2.2.2 we showed that choosing using MPE should be based on the position of the first peak after the noise in the MPE plot for an embedding dimension of . At this maximum, the normalized PE hits a maximum of approximately 1. From this methodology, we developed Algorithm 1 to determine the delay using the location of the first peak, while ignoring the noise region in Fig. 11.

       set ;
       start with a delay of ;
       start with initial normalized PE as ;
       while first peak not found do
             calculate normalized PE as ;
             if  then
                   Set outside of noise section flag as true ()
             end if
            if  and  then
                   if  then
                         First peak found;
                   end if
             end if
       end while
      return ;
Algorithm 1 Algorithm using MPE for .