On the Automatic Parameter Selection for Permutation Entropy

TL;DR

This paper introduces automated methods for selecting permutation entropy parameters, improving accuracy and efficiency in time series analysis across various systems by leveraging frequency-domain and information-theoretic approaches.

Contribution

It develops and compares new automatic parameter selection schemes for permutation entropy, extending existing methods and demonstrating their effectiveness across different types of systems.

Findings

Frequency approach accurately suggests delay for periodic and EEG data.

Mutual information with adaptive partitions best for chaotic systems.

False Nearest Neighbors and MPE reliably determine permutation dimension n.

Abstract

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 $\tau$. 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 $\tau$. We then compare our methods as well as current methods in the literature for obtaining both $\tau$ 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 $\tau$, 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.

Tables

Table 1. Table 1: A comparison between the calculated and suggested values for the delay parameter τ 𝜏 \tau for multiple MI approximation methods. The cells in bold highlight the methods that yielded the closest match to the suggested delay. The equal-sized partition method is described in Section  A.4.1 , Kraskov et al. methods 1 1 1 and 2 2 2 in Section  A.4.3 , and the adaptive partitioning approach in Section  A.4.2 .

	Mutual Information
System	Equal-sized Partitions	Kraskov et al. Method 1	Kraskov et al. Method 2	Adaptive Partitions	Suggested Delay $\mathrm{𝐭𝐚𝐮}$	Ref.
White Noise	1	3	3	1	1	[45]
Lorenz	13	9	9	9	10	[45]
Rossler	14	13	11	9	9	[51]
Bi-directional Rossler	16	14	14	15	15	[45]
Mackey-Glass	7	8	7	7	1 to 700	[45]
Sine Wave	4	17	13	1	15	[51]
Logistic Map	5	8	11	5	1 to 5	[45]
Henon Map	12	15	13	8	1 to 5	[45]
ECG	22	16	9	8	1 to 4	[45]
EEG	6	5	5	5	1 to 3	[45]

Table 2. Table 2: A comparison between the calculated and suggested values for the delay parameter τ 𝜏 \tau . The cells in bold show 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 ( 4 < n < 6 4 𝑛 6 4<n<6 ), AP under MI is an abbreviation for adaptive partitioning, and AC is an abbreviation for autocorrelation.

Catagory	System	Traditional Methods		Modified/Proposed Methods			Suggested Delay ($\tau$)	Ref.
		MI using AP	Spearman’s AC	Freq. App.	MPE	PAMI ($\mathrm{𝟒}\le 𝐧\le \mathrm{𝟔}$)
Noise	White Noise	1	1	1	1	1	1	[45]
Chaotic Differential Equation	Lorenz	9	15	6	17	5 to 9	10	[45]
	Rossler	9	12	7	19	6 to 10	9	[51]
	Bi-directional Rossler	15	12	7	20	6 to 10	15	[45]
	Mackey-Glass	7	5	3	8	2 to 4	1 to 700	[45]
Periodic	Sine Wave	1	10	21	16	5 to 8	15	[51]
Nonlinear Difference Eq.	Logistic Map	5	1	1	1	1	1 to 5	[45]
	Henon Map	8	1	1	1	1	1 to 5	[45]
Medical Data	ECG	8	21	2	13	1 to 2	1 to 4	[45]
	EEG	5	4	1	4	2 to 4	1 to 3	[45]

Table 3. Table 3: A comparison between the calculated and suggested values for the embedding dimension n 𝑛 n . The cells in bold show the methods that yielded the closest match to the suggested dimension.

Catagory	System	Traditional Methods		Modified Method	Suggested Dim. ($𝐧$)	Ref.
		FNN	SSA	MPE
Noise	White Noise	4	23	5	3 to 7	[45]
Chaotic Differential Equation	Lorenz	3	4	5	5 to 7	[45]
	Rossler	4	4	4	6	[51]
	Bi-directional Rossler	4	4	4	6 to 7	[45]
	Mackey-Glass	4	6	4	4 to 8	[45]
Periodic	Sine Wave	4	2	3	4	[51]
Nonlinear Difference Equation	Logistic Map	4	3	5	2 to 16	[45]
	Henon Map	4	2	5	3 to 10	[45]
Medical Data	ECG	7	8	5	3 to 7	[45]
	EEG	5	11	6	3 to 7	[45]

Equations

$$H(n)=-\sum p({\pi }_i)\log p({\pi }_i),$$

$$h_n=-\frac{1}{{\log }_{2}n!}\sum p({\pi }_i){\log }_{2}p({\pi }_i).$$

$$2f_{\max }<f_s<4f_{\max },$$

$$\tau =\frac{f_s}{\alpha f_{\max }},$$

$$\begin{array}{rl}LS& :\min \sum _i{r}_i^2,\\ LMS& :\min \left({\mathrm{median}}_i(r_i^2)\right),\end{array}$$

$$P_{|X|}(|X|)=\frac{2|X|}{E_w{\sigma }_x^2}{e}^{\frac{-|X{|}^2}{E_w{\sigma }_x^2}},$$

$$|X{|}_{\max }=\sqrt{\frac{E_w{\sigma }_x^2}{2}}.$$

$$C{P}_{|X|}(|X|)=1-e^{\frac{-|X{|}^2}{E_w{\sigma }_x^2}},$$

$$|X{|}_{\mathrm{cutoff}}=\sqrt{-E_w{\sigma }_x^2\ln (1-CP)}.$$

$$C=\frac{|X{|}_{\mathrm{cutoff}}}{|X{|}_{\max }}=\sqrt{-2\ln (1-CP)}.$$

$$d_2=\frac{1}{2}P=\frac{1}{f_s}{\tau }_r=\frac{1}{2f},$$

$$\tau =\frac{2}{\alpha }{\tau }_r,$$

$$\frac{dx}{dt}=\sigma (y-x),\frac{dy}{dt}=x(\rho -z)-y,\frac{dz}{dt}=xy-\beta z,$$

$$I(X;Y)=\sum _{x\in X}\sum _{y\in Y}p(x,y)\log \frac{p(x,y)}{p(x)p(y)},$$

$$I_p(\tau ,n)=H_{x(t,n)}+H_{x(t+\tau ,n)}-H_{x(t,n),x(t+\tau ,n)},$$

$$R_i=\sqrt{\frac{d_{m+1}^2(i,j)-d_m^2(i,j)}{d_m^2(i,j)}}.$$

$$D=\frac{g_L}{g_r}$$

$$h_n^{\prime }=\frac{-1}{n-1}H(n),$$

$$H(3)=-\frac{2}{5}\log \frac{2}{5}-\frac{2}{5}\log \frac{2}{5}-\frac{1}{5}\log \frac{1}{5}-\frac{1}{5}\log \frac{1}{5}=1.918bits.$$

$$m_{\mathrm{avg}}{\tau }_N\approx A_{\mathrm{noise}}{f}_s,$$

$${\rho }_{xy}=\frac{{\sum }_{i=1}^n(x_i-\bar{x})(y_i-\bar{y})}{\sqrt{{\sum }_{i=1}^n{(x_i-\bar{x})}^2}\sqrt{{\sum }_{i=1}^n{(y_i-\bar{y})}^2}}.$$

$$I(X;Y)=H(X)+H(Y)-H(X,Y),$$

$$I^{(1)}(X;Y)=\psi (k)-(\psi (n_x+1)+\psi (n_y+1)+\psi (N).$$

$$I^{(2)}(X;Y)=\psi (k)-1/k-[\psi (n_x)+\psi (n_y)]+\psi (N).$$

$$\frac{dx}{dt}=\sigma (y-x),\frac{dy}{dt}=x(\rho -z)-y,\frac{dz}{dt}=xy-\beta z.$$

$$\frac{dx}{dt}=-y-z,\frac{dy}{dt}=x+ay,\frac{dz}{dt}=b+z(x-c),$$

$$\begin{array}{rl}\frac{d{x}_1}{dt}& =-w_1{y}_1-z_1+k(x_2-x_1),\\ \frac{d{y}_1}{dt}& =w_1{x}_1+0.165y_1,\\ \frac{d{z}_1}{dt}& =0.2+z_1(x_1-10),\\ \frac{d{x}_2}{dt}& =-w_2{y}_2-z2+k(x_1-x_2),\\ \frac{d{y}_2}{dt}& =w_2{x}_2+0.165y_2,\\ \frac{d{z}_2}{dt}& =0.2+z_2(x_2-10),\end{array}$$

$$x(t)=-\gamma x(t)+\beta \frac{x(t-\tau )}{1+{x(t-\tau )}^n}$$

$$x(t)=sin(2\pi t)$$

$$x_{n+1}=r{x}_n(1-x_n),$$

Full text

Permutation Entropy (PE) is a cost effective tool for summarizing the complexity of a time series. It has been used in many applications including damage detection, disease forecasting, detection of dynamical changes, and financial volatility analysis. However, to successfully use PE, an accurate selection of two parameters is needed: the permutation dimension $n$ and embedding delay $\tau$. These parameters are often suggested by experts based on a heuristic or by a trial and error approach. Both of these methods can be time-consuming and lead to inaccurate results. In this work 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 Function (PAMI) and Multi-scale Permutation Entropy (MPE) for determining $\tau$. We then compare our methods as well as current methods in the literature for obtaining both $\tau$ 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 $\tau$, 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 a value of $n=5$ being suitable in most cases.

††preprint: AIP/123-QED

Permutation Entropy (PE) is a simple yet very effective tool for studying time series of dynamical systems. It provides an information statistic that measures the complexity of a time series through the probability of unique ordinal patterns found within the time series called permutations. These permutations are symbolic representations obtained by encoding consecutive subsets of the data of a certain length using their ordinal ranking. However, the success of PE depends on the selection of both the spacing (delay $\tau$) and size (dimension $n$) of these permutations. Despite the wide use of PE, it is often unclear how these parameters must be selected with the most common approach relying on trial and error. This can lead to inaccurate results, and it can prevent applying PE to large data sets in the absence of automatic parameter selection algorithms. In this work we investigate various methods for automatically selecting both $\tau$ and $n$. In addition to developing novel methods that facilitate the parameter selection, we also assess the accuracy of using classical time series tools for identifying permutation entropy parameters. The success for each of the investigated methods in computing $\tau$ and $n$ is determined based on a comparison with the corresponding values suggested by domain experts.

I Introduction

Permutation Entropy (PE) has its origins in information entropy, which is a tool to quantify the uncertainty in an information-based system. Information entropy was first introduced by Shannon shannon2001mathematical 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 $H_{s}(n)=-\sum{}p(x_{i})\log{p(x_{i})}$, where $x_{i}$ represents a data type, and $p(x_{i})$ is the probability of that data type. In recent years information entropy has been heavily applied to the time series of dynamic systems. Several new variations of information entropy have been proposed to better accommodate these applications, e.g. approximate entropy pincus1991approximate , sample entropy richman2000physiological , and PE bandt2002permutation . 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 effective for analyzing the dynamic state and complexity of a time series cao2004detecting ; myers2019persistent ; mccullough2015time ; Garland2018 ; Garland2018a ; Deng2015 ; Amigo2012 ; Bariviera2018 . PE is also noise robust for time series of sufficient length and relatively high signal-to-noise ratios, which is the ratio between useful signal and background noise. Alternatively, if the time series is relatively short or has a low signal-to-noise ratio, it is suggested to use a different entropy measurement such as coarse-grained entropies Porta2015 . 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 $\pi_{i}$. PE has two parameters: the permutation dimension $n$ and embedding delay $\tau$, which are used when selecting the permutation size and spacing, respectively. PE is sensitive to these parametersLi2012 ; STANIEK2007 ; riedl2013practical and there is no accurate selecting approach for all applications. 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) parameters suggested by experts for a specific application, (2) trial and error to find suitable parameters, or (3) methods developed for phase space reconstruction. We will now overview a simple example to better understand these parameters.

Bandt and Pompe bandt2002permutation defined PE according to

[TABLE]

where $p(\pi_{i})$ is the probability of a permutation $\pi_{i}$ and $H(n)$ is the permutation entropy of dimension $n$ with units of bits when the logarithm is of base 2. The permutation entropy parameters $\tau$ and $n$ are used when selecting the motif size, with $\tau$ determining the time difference between two consecutive points in a uniformly sub-sampled time series and $n$ as the permutation length or motif dimension. To form a permutation, begin with with an element $x_{i}$ of the series $X$. Using this element, the dimension $n$, and delay $\tau$, define the vector $v_{i}=[x_{i},x_{i+\tau},x_{i+2\tau},\ldots,x_{i+(n-1)\tau}]$. The corresponding permutation $\pi_{i}$ of this vector is determined using its ordinal pattern. For example, consider the third degree $n=3$ permutation shown in Fig. 1.

The permutation type, which categorizes the permutation, is found by first ordering the $n$ 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 $\pi_{i}=(1,0,2)$, which is one of $n!$ possible permutations for a dimension $n$, see Fig. 2 for the other possible permutations of $n=3$.

We can normalize PE using the maximum possible PE value, which occurs when all $n!$ possible permutations are equiprobable according to $p(\pi_{1})=p(\pi_{2})=\ldots=p(\pi_{n!})=\frac{1}{n!}$. The resulting normalized PE is

[TABLE]

A toy example demonstrating the calculation of $h_{n}$ is provided in the appendix.

Many domain scientists who apply PE make general suggestions for $n$ and $\tau$ zhang2018analysis ; frank2006permutation , which can be impractical for some applications. As an example, Popov et al. popov2013permutation showed the influence of the sampling frequency on the proper selection of $\tau$. As for the dimension $n$, there are general suggestions riedl2013practical on how to choose its value based on the vast majority of applications having an appropriate permutation dimension in the range $3<n<8$. Additionally, Bandt and Pompe bandt2002permutation suggest that $N\gg n$, where $N$ is the length of the time series. However, these general outlines for the selection of $n$ (and $\tau$) do not allow for an application specific suggestions.

If we assume that suitable PE parameters correspond to optimal phase space reconstruction parameters, then a common approach for selecting $\tau$ and $n$ is to implement one of the existing methods for estimating the optimal Takens’ embedding takens1981detecting parameters. Hence, some of the common methods for determining $\tau$ include the mutual information function approach fraser1986independent , the first folding time of the autocorrelation function grassberger1983measuring ; Box2015 , and phase space methods buzug1992optimal . Additionally, some common phase space reconstruction methods for determining $n$ include box-counting block1990efficient , correlation exponent method grassberger1983measuring , and false nearest neighbors kennel1992determining . 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 also lead to inaccurate $n$ or $\tau$ values. In spite of these differences, permutations can be viewed as symbolic representation of regions in the phase space through a binning process. Permutations partition the phase space based on the ordinal rankings of the embedded vectors. This relationship between phase space and permutations opens up the potential for some of the classic phase space reconstruction methods for selecting both $n$ and $\tau$ to be a plausible solution for selecting the same parameters for PE.

Even with the possibility that phase space reconstruction methods for selecting $\tau$ and $n$ may work for choosing synonymous parameters of PE, there are a few practical issues that preclude using parameters from time series reconstruction for PE. One issue stems from many of the methods (e.g. false nearest neighbors and mutual information) still requiring some degree of user input through either a parameter setting or user interpretation of the results. This introduces issues for practitioners working with numerous data sets or those without enough expertise in the subject area to interpret the results. Another issue that arises in practice is that the algorithmic implementation of existing time series analysis tools is nontrivial. This hinders these tools from being autonomously applied to large datasets. For example, the first minimum of the MI function is often used to determine $\tau$. However in practice there are limitations to using mutual information to analyze data without the operator intervention to sift through the minima and choose the first ’prominent’ one. This is due to possibility that the mutual information function can have small kinks that can be erroneously picked up as the first minimum. Figure 3a shows this situation, where the first minimum of the mutual information function for a periodic Lorenz system is actually an artifact and the actual delay should be at the prominent minimum with $\tau=11$. Further, the mutual information function approach may also fail if the mutual information is monotonic. This is a possibility since there is no guarantee that minima exist for mutual information Babaie-Zadeh2005 . An example of this mode of failure is shown in Fig. 3b, which was generated using EEG data Andrzejak2001 from a patient during a seizure.

A mode of failure for the autocorrelation method can occur when the time series is non-linear or has a moving average. In this case, the autocorrelation function may reach the folding time at an unreasonably large value for $\tau$. As an example, Fig. 3c shows the autocorrelation not reaching the folding time of $\rho=1/e$ until a delay of $\tau=283$ for electrocardiogram data provided by the MIT-BIH Arrhythmia Database Moody2001 . The last mode of failure concerns choosing the permutation dimension $n$ to be equal to the embedding dimension optimized using delay embedding from time series analysis. This can lead to an overly large embedding dimension Chung2019 ($n\gg 8$), which would make the calculation of PE impractical because the number of possible permutations $n!$ would become too large.

All of these possible modes of failure can make using classical phase space methods for selecting $\tau$ and $n$ unreliable thus necessitating new tools or modifications to make selecting $\tau$ and $n$ for PE more robust and less user-dependent.

These shortcomings lead us to the problem that we address in this paper: Given a sufficiently sampled/oversampled and noisy time series $X=\{x_{t}\}_{\mathcal{R}^{+}}$, how can we reliably and systematically define appropriate dimension $n$ and time delay $\tau$ values for computing the corresponding PE?

Our first contribution towards answering this question is detailed in Section II, which addresses the automatic selection of the time delay $\tau$. In Section II.1 we combine the Least Median of Squares (LMS) approach for outliers detection with Fourier transformation theorem to derive a formula for the maximum significant frequency in the Fourier spectrum, with the assumption that $X$ is contaminated by Gaussian measurement noise. This formula allows obtaining a cutoff value where the only input, besides the time series, 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 $\tau$ value.

The second contribution is through an approach that we develop in Section II.2, which uses Multi-scale Permutation Entropy (MPE) for finding $\tau$. 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 $\tau$ for a chaotic time series by using the first maxima in the MPE as it satisfies the Nyquist’s sampling theorem.

Our third contribution to the automatic selection of $\tau$ is through the analysis of Permutation Auto-Mutual Information Liang2013 (PAMI). PAMI is an existing method for measuring the mutual information of permutations. However, we tailor this method to specifically select $\tau$ for PE.

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 $n$. We compare dimension $n$ values computed from False Nearest Neighbors (FNN—Section III.1), Singular Spectrum Analysis (SSA—Section III.2), and MPE (Section II.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 riedl2013practical , which requires finding $\tau$ first. We made this process automatic through the selection of $\tau$ from our second contribution.

This paper is organized as follows. We first go into detail on some existing methods for selecting both $\tau$ and $n$. Specifically, in Section II we provide a detailed explanation for selecting $\tau$ using existing, automatic methods such as autocorrelation in Section II.3 and Mutual Information (MI) in Section II.4. Additionally, we modify and develop/tailor methods to automatically select $\tau$. These methods include a frequency approach in Section II.1, MPE in Section II.2, and PAMI in Section II.5. In Section III we expand on the process for selecting $n$ using False Nearest Neighbors (FNN) in Section III.1 and Singular Spectrum Analysis in Section III.2. In Section III.3, we explain our algorithm for automatically selecting $n$ using MPE. After introducing each method, in Section IV we contrast all of these methods and make conclusions on their viability by comparing the resulting parameters to those suggested by PE experts. An overview of the methods that will be investigated for automatically calculating both $\tau$ and $n$ are shown in Fig. 4. All the functions used and developed in this work are available in Python through GitHub Myers2020 .

II Methods for Embedding Delay Parameter Selection

The delay embedding parameter $\tau$ is used to uniformly subsample the original time series. To elaborate, consider the time series $X=\{x_{i}\mid i\in\mathbb{N}\}$. By applying the delay $\tau\in\mathbb{N}$, a new subsampled series is defined as $X(\tau)=[x_{0},x_{\tau},x_{2\tau},\ldots]$. In order to obtain a stable and automatic method for estimating an optimal value for $\tau$ we investigate: a novel frequency-based analysis that we describe in Section II.1, Multi-scale Permutation Entropy (MPE) (Section II.2), autocorrelation (Section II.3), and Mutual Information function (MI) (Section II.4). We recognize, but do not investigate, some other methods for finding $\tau$ such as diffusion maps Berry2013 and phase space expansion buzug1992optimal .

II.1 Frequency Approach for Embedding Delay

In this section we develop a method for finding the noise floor in the Fourier spectrum using Least Median of Squares (LMS) massart1986least . 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 melosik20160 , which agrees with Nyquist sampling theorem landau1967sampling for choosing a suitable sampling frequency $f_{s}$ as

[TABLE]

where $f_{\rm max}$ is the maximum significant frequency in the signal. Melosik and Marszalek melosik20160 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. However, the automatic identification of $f_{\rm max}$ from an oversampled signal is not trivial. Melosik and Marszalek melosik20160 selected a maximum significant frequency by inspecting the normalized Fourier spectrum and using a threshold cutoff of approximately $0.01$ 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 $f_{\rm max}$. Further, attempting to algorithmically adopt the approach suggested by Melosik and Marszalek melosik20160 resulted in large errors especially in the presence of a low signal to noise ratio. This motivated the search for an automatic and data-driven approach for identifying the noise floor which could then be used to find the maximum significant frequency. To do this we develop a method based on $1$-D least median of squares applied to the Fourier spectrum. The assumptions inherent to our method are

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 time series is undersampled. 2. 2.

The Fourier transform of the time series needs to have less than 50$\%$ of the points with significant amplitudes. This requirement stems from the limitations of the least median of squares regression. 3. 3.

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

We find suitable cutoffs for obtaining $f_{\rm max}$ of the signal by using the noise floor determined from the $1$-D least median of squares, and compute a suitable embedding delay according to

[TABLE]

where we set $\alpha=2$, thus agreeing with the range in Eq. (3) and the Nyquist sampling criterion.

Figure 5 summarizes the frequency approach for $\tau$ with the use of our 1-D LMS method for finding a noise floor in the Fourier spectrum. This process begins with computing the Fourier spectrum of the signal, which is followed by fitting an 0-D LMS regression line to the noise in the Fourier spectrum. This provides statistical information about the Probability Distribution Function (PDF) of the noise level. The PDF is used to determine the Cumulative Distribution Function (CDF), which we use determine a meaningful noise cutoff in the Fourier spectrum. However, it is assumed that the noise is approximately GWN for this method to hold statistical significance. This cutoff is used to separate the highest significant frequency in the Fourier spectrum $f_{\rm max}$, which is used to find a suitable embedding delay $\tau$ 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.

Least Median of Squares:

LMS massart1986least is a robust regression technique used when up to 50$\%$ of the data is corrupted by outliers. Outliers will be considered as anything other than noise in the fourier spectrum for our application. 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 is highlighted as

[TABLE]

where $r$ is the residual. Similar to the $i$ subscript in $\sum_{i}$, the $i$ in ${\rm median}_{i}$ signifies that the median is of all residuals. Figure 6 shows an example application of the linear LMS regression.

Specifically, this figure shows 110 data points drawn from the line $y=x+1$ with added GWN of zero mean and 0.1 standard deviation. The data is corrupted with 90 outliers centered around $(3,2)$ with a normal distribution of 1.0 along $x$ and 0.6 along $y$. Figure 6 shows that the linear regression results closely match the actual trend line with the fitted line being $y=0.998x+1.012$ in comparison to the actual $y=x+1$.

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 white noise. Additionally, this section derives the location of the theoretical maximum of the PDF. The FT distribution of GWN richards2013discrete is described as

[TABLE]

where $|X|$ is the magnitude of the FT of GWN, $P_{|X|}$ is the probability density function of $|X|$, $\sigma_{x}$ is the standard deviation of the GWN, and $E_{w}$ is the window energy or number of discrete transforms taken during the FT. By setting the first derivative of $P_{|X|}$ with respect to $|X|$ equal to zero, the theoretical maximum of the PDF is

[TABLE]

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 as papoulis2002probability

[TABLE]

where $CP_{|X|}$ is the cumulative probability of $|X|$.

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 0-D fit of the Fast Fourier Transform (FFT) of the signal, which results in an approximate value of $|X|_{\rm max}$, which is $|X|$ at the maximum of $P_{|X|}$. Using $|X|_{\rm max}$ from the LMS fit, we then find the standard deviation of the distribution $\sigma_{x}$ 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 $|X|_{\rm max}$. 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. 8A) results in the distribution shown in Fig. 8B. 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 $|X|_{\rm max}$. Additionally, the theoretical shape of the PDF in Fig. 8B is shown to be very similar to the actual distribution.

Next, our approach utilizes Eq. (8) and $\sigma_{x}$ derived from Eq. (7) for finding the cutoff value $|X|_{\rm cutoff}$. The $|X|_{\rm cutoff}$ for a desired cumulative probability $CP$ is found by solving Eq. (8) for $|X|$ as

[TABLE]

In order to make $|X|_{\rm cutoff}$ robust to normalization and scaling of the FFT we define the ratio $C$ between the suggested cutoff from Eq. (9) and the maximum of the PDF from Eq. (7) as

[TABLE]

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 $CP=99\%$ 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 suggest a cutoff ratio $C=6$ to be used for signals with less than $10^{4}$ data points. This yields an expected probability of $\approx 10^{-8}\%$ for a point in the FFT of the GWN attaining a magnitude greater than $|X|_{\rm cutoff}$. Alternatively, Eq. (10) can be used to calculate a different value of C based on the desired probability and length of the signal.

II.2 Multi-scale Permutation Entropy for Selecting 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 for analyzing physiological time series costa2002multiscale . Zunino et al. zunino2010permutation showed how the first maxima in the MPE plot arises when $\tau$ matches the characteristic time delay $\tau_{r}$. Furthermore, the periodicity can be captured by the first dip in the MPE plot as shown in Fig. 9 at the location $\rm d_{2}$ when the delay $\tau$ matches the characteristic time delay $\tau_{r}$.

Figure 9 shows embedding delays $d_{0}$, $d_{1}$, and $d_{2}$ calculated as $d=\frac{\tau}{f_{s}}$ as well as their corresponding locations on a normalized MPE plot. This toy MPE plot shows that the normalized MPE reaches its first maximum when the delay is roughly $d_{1}$, which corresponds to approximately an even distribution of permutations. A second observation, as mentioned previously, is that at $d_{2}$ (or the first dip in the MPE plot) there is a resonance or aliasing effect caused by $\tau\approx\tau_{r}$, which can be used to determine the period of the time series. This is based on the embedding delay size at $d_{2}$ causing the embedding vector size $V=d(n-1)$ to be approximately half the of the periodicity P, which can be expressed as

[TABLE]

where $P$ is the main period of oscillation, $f$ is the main frequency of the time series corresponding to $P$, and $f_{s}$ 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 reduces PE through more regularity in the permutation distribution.

We use the criteria of Melosik and Marszalek melosik20160 to determine a suitable delay from the location of the first dip at $d_{2}$. 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 $\tau$. However, for MPE, we substitute $f_{s}$ and $f_{\rm max}$ in Eq. (3) with $f_{s}=2f\tau_{r}$ from Eq. (11) and $f_{\rm max}=f$. These substitutions allow Eq. (4) to reduce to

[TABLE]

where $\alpha\in[2,4]$. These simplifications show that $\tau$ is only dependent on the delay which causes resonance $\tau_{r}$ when applying MPE. However, for a chaotic time series, the dip at $\tau_{r}$ may not be present due to non-linear trends. To address this issue, we will first investigate the three dominant regions of the MPE plot, which will also be located for a chaotic time series example. We will then propose a new, automatic method for selecting $\tau$ that agrees with the frequency criteria stated in Eq. (12). Additionally, in Section A.2.2 of the appendix we investigate the robustness of the method to noise and in Section A.5 of the appendix we provide the algorithm (Algorithm 1) for finding $\tau$ using MPE.

MPE Regions

Riedl et al. riedl2013practical showed that the MPE plot can be separated into three distinct regions as described below and shown in Fig. 10. 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 $\tau$ causes the motif distribution corresponding to the permutation entropy to be heavily weighted on just increasing or decreasing motifs (motifs (0,1,2) and (2,1,0) for $n=3$ from Fig. 2). This effect was coined as the “Redundancy Effect" by De Micco et al. de2012sampling , which means sufficiently low values of $\tau$ result in redundant motifs. However, as $\tau$ 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 shows a slight dip to the first minima. This reduction in permutation entropy is caused by the aliasing or resonance from the value of $d$ 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 has possible additional minima and maxima from additional alignment of the embedding vector $d$ with multiples of the main period. This region was referred to as the “Irrelevant Region" by De Micco et al. de2012sampling due to effectively large values of $\tau$ forcing the delayed sampling frequency to fall below the Nyquist sampling rate as described by the lower bound in Eq. (3).

MPE Example with Chaotic Time Series

In Sections II.2 and II.2, 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 $\tau$. In this section we further show the applicability of this approach to chaotic signals using the $x$-coordinate of the Lorenz System as an example. We simulate the Lorenz equations

[TABLE]

with a sampling rate of $100$ Hz and using the parameters $\rho=28.0$, $\sigma=10.0$, and $\beta=8.0/3.0$. This system was solved for $100$ seconds and only the last $15$ seconds from the time series are used. Figure 11 shows the result of applying MPE to the simulated Lorenz system.

Figure 11 shows similarities to Fig. 10 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 $\tau_{r}$. We suggest using the first maxima to find $\tau$ because this delay is likely to fall within the region described by Eq. (12).

II.3 Autocorrelation for Embedding Delay

Autocorrelation is a traditional method for selecting $\tau$ for phase space reconstruction by using the correlation coefficient between the time series and its $\tau$-lagged version. This method was first introduced by Box et al. Box2015 . Typically, the autocorrelation function is computed as a function of $\tau$ and, as a rule of thumb, a suitable delay $\tau$ is found when the correlation between $x(t)$ and $x(t+\tau)$ reaches the first folding time, i.e., when $\rho\leq 1/e$ Kantz2004 . The two prominent correlation techniques that are commonly used when implementing an autocorrelation-based approach for finding $\tau$ are Pearson Correlation (see Section A.3.1 of appendix) and Spearman’s Correlation (see Section A.3.2 of appendix). Additionally, an example demonstrating how to calculate $\tau$ using autocorrelation and the difference between the two correlation methods is provided in Section A.3.3 of the appendix.

II.4 Mutual Information for Embedding Delay

Mutual information (MI) can be used to select the embedding delay $\tau$ 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. shannon1951mathematical as

[TABLE]

where $X$ and $Y$ are the two sequences, $p(x)$ and $p(y)$ are the probability of the element $x$ and $y$ separately, and $p(x,y)$ is the joint probability of $x$ and $y$. Fraser and Swinney fraser1986independent showed that for a chaotic time series the MI between the original sequence $x(t)$ and and delayed version $x(t+\tau)$ will decrease as $\tau$ increases until reaching a first minimum. At this minima, the delay $\tau$ 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 used for the selection of the PE parameter $\tau$. 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 $\tau$ using this method. All MI methods can be applied to either ranked or unranked data. We investigate four methods for estimating $\tau$ for PE using MI. These methods include MI with equal-sized partitions, adaptive partitions, and two permutation-based MI estimation methods. For details on these methods please reference the appendix in Section A.4.

To determine the optimal MI approximation method for selecting $\tau$ for PE, Fig. 12 shows a comparison between the $\tau$ values computed from each of the MI methods and the corresponding values suggested by experts. The table shows that the adaptive partitioning method of Section A.4.2 results in an accurate selection of $\tau$ for the majority of systems. We will use the adaptive partitioning estimation method when making comparisons to other methods. For the exact values of $\tau$ from each of the MI methods please reference Table 1 in the appendix.

II.5 Permutation Auto-mutual Information for Selecting Delay

As shown in Section II.4, Mutual information (MI) is a useful method for selecting $\tau$ for phase space reconstruction. However, it does not account for the permutation distribution when selecting $\tau$, which can lead to inaccuracies in computing the PE. To circumvent this issue, we develop a new method for selecting $\tau$ using Permutation Auto-Mutual Information (PAMI) Liang2013 , which was developed to detect dynamic changes in brain activity. We are tailoring PAMI for its application in the selection of the permutation entropy parameter $\tau$ for the first time. This is done by measuring the joint probability between the original permutations formed when a delay of $\tau=1$ is used and to the permutations when $\tau$ is incremented. PAMI is defined as

[TABLE]

where $H$ is the permutation entropy described in Eq. (1). We suggest an optimal delay $\tau$ for a given dimension $n$ when PAMI is at a minimum. This delay corresponds to minimum shared information between the original permutations with $\tau=1$ and its time lagged permutations. By applying this method for the simple sinusoidal function described in Section A.5.5, we can form Fig. 13 with $n\in[2,5]$ and $\tau\in[1,50]$.

As shown, the window size is approximately independent of the dimension $n$, with an optimal window $\tau(n-1)\approx 25$ 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. Hence, we implement the common choice for $n$ ranging from $4\leq n\leq 6$ for PE riedl2013practical . To reduce the computational demand, we suggest using permutation dimensions $n=2$ to find an optimal window size. In addition to the reduced computational demand of using $n=2$, we found that $I_{p}(n=2)\approx 0$ at the first minima. This also helps making this first minima even more simple.

III Methods for Motif Dimension

The second parameter for permutation entropy that needs to be automatically identified is the embedding dimension $n$. The methods for determining $n$ fall into one of two categories: (1) independently determining $n$ and $\tau$, and (2) simultaneously determining $n$ and $\tau$ based on the width of the embedding window. For the first category, we investigate using the method of False Nearest Neighbors (FNN) kennel1992determining in Section III.1, and Singular Spectrum Analysis (SSA)broomhead1986extracting in Section III.2. For the second category, we contribute to the selection of $n$ by developing an automatic method using MPE from Section III.3. This method combines the results for finding $\tau$ through MPE in Section II.2 with the work of Riedl et al. riedl2013practical . We acknowledge that our work does not include other commonly used methods for independently calculating $n$ such as box-counting clark1990estimating , largest Lyapunov exponent wolf1985determining , and Kolmogorov–Sinai entropy pesin1977characteristic .

III.1 False Nearest Neighbors for Embedding Dimension

False Nearest Neighbors (FNN) is one of the most commonly used methods for geometrically determining the minimum embedding dimension $n$ for state space reconstruction  kennel1992determining . For this method the time series is repeatedly embedded into a sequence of $m$-dimensional Euclidean spaces for a range of increasing values of $m$. The idea is that when the minimum embedding dimension $m$ is reached or $m\geq n$, the distance between neighboring points does not significantly change as we keep increasing $m$. In other words, the Euclidean distance $d_{m}(i,j)$ between the point $\mathbf{P}_{i}\in\mathbb{R}^{m}$ and its nearest neighbor $\mathbf{P}_{j}\in\mathbb{R}^{m}$ minimally changes when the embedding dimension increases to $m+1$. If the dimension $m$ is not sufficiently high, then the points are false neighbors if their pairwise distance significantly increases when incrementing $m$. This ratio of change in the distance between nearest neighbors embedded in $\mathbb{R}^{m}$ and $\mathbb{R}^{m+1}$ is quantified using the ratio of false nearest neighbors

[TABLE]

$R_{i}$ is compared to the tolerance threshold $R_{\rm tol}$ to distinguish false neighbors when $R_{i}>R_{\rm tol}$. In this paper, we select $R_{\rm tol}=15$ as used by Kennel et al. kennel1992determining . By applying this threshold over all points, we can find the number of false neighbors as a percent FNN $P_{\rm FNN}$. If there is no noise in the system, $P_{\rm FNN}$ should reach zero when a sufficient dimension is reached. However, with additive noise present, $P_{\rm FNN}$ may never reach zero. Thus, it is commonly suggested to use a percent FNN cutoff for finding a sufficient dimension $n$. We use the typically chosen cutoff $P_{\rm FNN}<10\%$, which is suitable for most applications when moderate noise is present.

III.2 Singular Spectrum Analysis for Embedding Dimension

The singular spectrum analysis method was first introduced as a tool to find trends and prominent periods in a time seriesbroomhead1986extracting . Leles et al. leles2018improving 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 $L$ 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 $D$

[TABLE]

of the sum of the $L$th diagonal entries $g_{L}$ to the sum of the total diagonal entries $g_{r}$. When $D$ exceeds $0.9$, we consider the dimension to be high enough and set $n=L$, which can then be used as the embedding dimension for permutation entropy.

III.3 Multi-scale Permutation Entropy for Permutation Dimension

Riedl et al. riedl2013practical showed how MPE can be used to determine an embedding dimension $n$. This method requires the embedding delay $\tau$ to be set to the length of the main period of the signal as shown in Section II.2. The theory behind the method is based on normalizing the MPE according to

[TABLE]

where $h^{\prime}_{n}$ is the PE normalized using the embedding dimension, and $H_{n}$ is the PE calculated from Eq. (1). Riedl et al. riedl2013practical determine the embedding dimension by incrementing $n$ to find the largest corresponding normalized PE $h^{\prime}_{n}$ with an embedding delay $\tau$ heuristically determined from the main period length. They concluded that the $h^{\prime}_{n}$ with the highest entropy accurately accounts for the needed complexity of the time series, and therefore suggests a suitable embedding dimension. Rield et al. riedl2013practical 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 $\tau$.

To make the process automatic, we introduce algorithm 1 based on Section II.2 to automatically select the correct $\tau$, which we then use in conjunction with Eq. (18) to find $n$ corresponding to the maximum $h^{\prime}_{n}$. Additionally, we suggest scaling $n$ from $3$ to $8$ as we have not yet found a system requiring $n>8$ using this method.

IV Results and Discussion

To make conclusions about the described methods for determining $\tau$ and $n$, we made comparisons to values suggested by experts. The majority of the suggested parameters are taken from the work of Riedl et al. riedl2013practical , while parameters for the Rossler system and sine wave are from Tao et al. tao2018permutation . Figures 14 and 15 show the calculated and suggested values for $\tau$ and $n$, respectively. For the exact values of $\tau$ and $n$ from each of the parameter estimation methods please reference Tables 2 and 3 in the appendix, respectively. Additionally, script for reproducing the results found in this paper are provided through the Mendeley.

Embedding Delay

Figure 14 shows the automatically computed $\tau$ in comparison to the expert-identified values for a variety of systems. These systems fall within several categories including the following: noise, chaotic differential equations, periodic systems, nonlinear difference equations, and medical data. The methods presented in Fig. 14 include PAMI from Section II.5, MI calculated using adaptive partitioning from Section A.4.2, Spearman’s Autocorrelation from Section II.3, MPE from Section II.2, and the frequency approach from Section II.1. For the noise category we only investigated Gaussian white noise, 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. However, as addressed in Section I, there are possible modes of failure for MI. To validate that MI is accurately selecting a value for $\tau$, we recommend also calculating $\tau$ using the frequency approach. 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. Therefor, we suggest using both of these methods to calculate $\tau$ for periodic systems. Additionally, we do not suggest the use of MI for periodic systems as it can have early false minima resulting in inaccurate delay selection. 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 $\tau$. As a generalization of the results found, we suggest the use of MI with adaptive partitioning when selecting $\tau$ for chaotic differential equations. For periodic systems, nonlinear difference equations, and ECG/EEG data we suggest the use of the frequency approach that we developed in this paper. However, when applying the frequency approach to quasiperiodic time series with multiple harmonics of decreasing amplitude, the method may fail due to the delay being selected based on an insignificant high frequency. The use of either Spearman’s autocorrelation or MPE may be more suitable under this condition. In general, multiple methods should be used for each system to validate that an accurate delay is selected due to the possible modes of failure of each method. Specifically, The frequency approach may fail if the noise does not have a Gaussian distribution, MI can fail if a false minima occurs or the relationship is monotonic, and autocorrelation can fail if the time series being analyzed does not oscillate about a fixed value.

Embedding Dimension

Figure 15 shows the automatically computed parameter $n$ in comparison to the expert-identified values. 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 $n$. 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. However, FNN may fail if the effects of noise are not correctly accounted for, which can lead to overly large embedding dimensions. These results also show that the dimension $n=6$ works well for almost all applications.

V Conclusion

In this paper we demonstrated various methods for automatically determining the PE parameters $\tau$ and $n$ 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 $\tau$ 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 $n$ include false nearest neighbors, singular spectrum analysis, and multiscale permutation entropy. Several of these methods for calculating $\tau$ or $n$ do have suggested parameters to be set by hand. This leaves some methods as not completely automatic. However, the methods of MI, autocorrelation, MPE, and PAMI do not have any parameters set, which reduces the user influence on parameter selection and improves the automatic selection. Additionally, the parameters that are suggested are default parameters that work for the majority of applications.

Our first contribution was developing a new frequency approach analysis and extending two existing methods, PAMI, and MPE, to automatically determine $\tau$. 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 using PAMI to select $\tau$, we showed how the first minimum in the PAMI can be used to find an optimal embedding window $n(\tau+1)$, where $\tau$ was then selected using the range $4\leq n\leq 6$. 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 $\tau$. 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 $\tau$ and $n$ using MPE. We also collected and compared some of the most popular methods for obtaining $n$ 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 $\tau$ and $n$. We found that SSA did not provide suitable values for $n$. However, both FNN and MPE provided accurate values for $n$ for most of the systems. We conclude that, for the majority applications, a permutation dimension $n=5$ is suitable. For determining $\tau$, we showed that our frequency approach provided accurate suggestions for $\tau$ 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.

Acknowledgment

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

Appendix A Appendix

A.1 PE Calculation Example

To demonstrate how PE is calculated, consider the sequence $X=[4,7,9,10,6,11,3,2]$ with PE parameters $n=3$ and $\tau=1$. The left side of Fig 16 shows how the sequence can be broken down into the following permutations: two $(0,1,2)$, one $(1,0,2)$, two $(1,2,0)$, and one $(2,1,0)$ 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. 16 shows the abundance of each permutation.

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

[TABLE]

This summarizes a simple example of how to apply permutation entropy to time series data.

A.2 MPE Algorithm and Effects of Noise

A.2.1 MPE delay Algorithm

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

A.2.2 Effects of Noise

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. 10, while only significantly affecting region A as shown in Fig 17. Furthermore, depending on the signal to noise ratio, there will only be an effect at the beginning of region A. Figure 17 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

[TABLE]

where $m_{\rm avg}$ is the average of the absolute value of the slope and $A_{\rm noise}$ is approximately the maximum amplitude of the noise, $\tau_{N}$ is the value of $\tau$ great enough to surpass the noise amplitude. We derived this condition from the need for, on average, $|f(t)-f(t+\tau)|>A_{\rm noise}$. This shows that MPE is robust to noise as long as the noise amplitude does not exceed the amplitude of the signal.

A.3 Autocorrelation Methods and Example

A.3.1 Pearson Correlation

The Pearson correlation coefficient $\rho_{xy}\in[-1,1]$ measures the linear correlation of two time series $x$ and $y$. Using these two data sets the correlation coefficient is calculated as

[TABLE]

The possible values of $\rho_{xy}$ represent the relationship between the two data sets, where $\rho_{xy}=1$ represents a perfect positive linear correlation, $\rho_{xy}=0$ represents no linear correlation, while $\rho_{xy}=-1$ 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.

A.3.2 Spearman’s Correlation

Spearman’s correlation is also calculated using Eq. (20) with the substitution of $x$ and $y$ 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. 18 shows two sequences $x$ and $y$ calculated from $y=x^{4}$ with $x\in[0,10]$. Using this example, the Pearson correlation is calculated as $\rho\approx 0.86$, while Spearman’s ranked correlation yields $\rho=1.0$. This result demonstrates how Spearman’s correlation coefficient accurately detects the non-linear, monotonic correlation between $x$ and $y$ whereas Pearson correlation may miss it.

A.3.3 Autocorrelation Example

We can use the concept of correlation to select a delay $\tau$ by calculating the correlation coefficient using Eq. (20) between a time series and its $\tau$-lagged version. As an example, take the time series $x(t)=\sin(2\pi t)$, with $t\in[0,5]$ having a sampling frequency of $100$ Hz. This results in a suggested delay $\tau=20$ at the first folding time using both Spearman’s and Pearson correlation. In section IV we will implement Spearman’s version of autocorrelation to account for the possibility of non-linear correlations.

A.4 MI methods

A.4.1 MI using Equal-sized Partitions

For the calculation of MI, the joint and independent probabilities of the original $x(t)$ and time lagged $x(t+\tau)$ time series are needed. However, since $x$ 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 $k$. We investigated various methods for estimating an appropriate number of bins using the length of the time series $N$. These methods include the common square-root choice $k=\lceil\sqrt{N}\rceil$, Sturge’s formula [49] $k=\lceil{\log_{2}({N})}\rceil+1$, and Rice Rule [28] $k=\lceil 2N^{1/3}\rceil$. After comparing each method using a variety of examples, we found that the use of Sturge’s formula provided the best results for selecting $\tau$ for PE using MI.

A.4.2 MI using Adaptive Partitions

Darbellay and Vajda [16] introduced a multistep, adaptive partitioning scheme to select appropriate binning sizes in the observation space formed by the plane $x(t)$ and $x(t+\tau)$. Their method is often considered state-of-the-art for estimating the mutual information function [26]. In this approach, the bins are recursively created where in the first function call, the space of the signal and its $\tau$-lagged version is divided into an equal number of $2$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. (14).

A.4.3 Kraskov MI

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

[TABLE]

where $H$ is the Shannon entropy. Next, an approximation of $H(X)$ with digamma functions is done, but the probability density of $X$ and $Y$ still needs to be estimated. To do this, adaptive partitions using the $k$-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 $k=1$ nearest neighbor $j$ to form square bins as shown in Fig. 19-a, and the second method in Fig. 19-b uses rectangular partitions using the horizontal and vertical distances to the $k=1$ nearest neighbor $j$.

To continue with the example shown in Fig. 19, the density probability is estimated using the strips formed from these bins. To highlight the difference, Fig. 19-a shows a horizontal strip of width $\epsilon(i)$ encapsulating $n_{x}(i)=2$ points (strip does not include the point $i$), while in Fig. 19-b only $n_{x}(i)=1$ point is enclosed. Using these probability density approximations and the digamma function $\psi$, MI between $X$ and $Y$ can be estimated. Using the partitioning method shown in Fig. 19-a the MI is estimated as

[TABLE]

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

[TABLE]

For the results shown in Section IV, we use the $k=3$ nearest neighbor to generate the partitions.

A.5 Dynamic System Models

A.5.1 Lorenz System

The Lorenz system used is defined as

[TABLE]

The Lorenz system had a sampling rate of 100 Hz with parameters $\sigma=10.0$, $\beta=8.0/3.0$, and $\rho=95$. This system was solved for 100 seconds and the last 24 seconds were used.

A.5.2 Rössler System

The Rössler system used was defined as

[TABLE]

with parameters of $a=0.1$, $b=0.1$, $c=14$, 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.5.3 Bi-Directional Coupled Rössler System

The Bi-directional Rössler system is defined as

[TABLE]

with $w_{1}=0.99$, $w_{2}=0.95$, and $k=0.05$. 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.5.4 Mackey-Glass Delayed Differential Equation

The Mackey-Glass Delayed Differential Equation is defined as

[TABLE]

with $\tau=2$, $\beta=2$, $\gamma=1$, and $n=9.65$. 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.5.5 Periodic Sinusoidal Function

The sinusoidal function is defined as

[TABLE]

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

A.5.6 EEG Data

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

A.5.7 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 [35]. We used data points 3000 to 4500 during normal sinus rhythm.

A.5.8 Logistic Map

The logistic map was generated as

[TABLE]

with $x_{0}=0.5$ and $r=3.95$. Equation 29 was solved for the first 500 data points.

A.5.9 Hénon Map

The Hénon map was solved as

[TABLE]

where $b=0.3$, $x_{0}=0.1$, $y_{0}=0.3$, and $a=1.4$. This system was solved for the first 500 data points of the x-solution.

A.6 Tabulated PE parameters