Confidence Intervals and Hypothesis Testing for the Permutation Entropy with an application to Epilepsy

TL;DR

This paper introduces a bootstrap-based statistical method to accurately estimate the variability of Permutation Entropy in time series, with applications to detecting epileptic activity in EEG signals.

Contribution

It proposes a novel parametric bootstrap approach for assessing the accuracy of Permutation Entropy, addressing a gap in statistical analysis of this complexity measure.

Findings

The method performs well on stochastic process simulations.

It provides reliable confidence intervals for Permutation Entropy estimates.

Application to EEG data distinguishes normal from pre-ictal states.

Abstract

In nonlinear dynamics, and to a lesser extent in other fields, a widely used measure of complexity is the Permutation Entropy. But there is still no known method to determine the accuracy of this measure. There has been little research on the statistical properties of this quantity that characterize time series. The literature describes some resampling methods of quantities used in nonlinear dynamics - as the largest Lyapunov exponent - but all of these seems to fail. In this contribution we propose a parametric bootstrap methodology using a symbolic representation of the time series in order to obtain the distribution of the Permutation Entropy estimator. We perform several time series simulations given by well known stochastic processes: the 1=f? noise family, and show in each case that the proposed accuracy measure is as efficient as the one obtained by the frequentist approach of repeating the experiment. The complexity of brain electrical activity, measured by the Permutation Entropy, has been extensively used in epilepsy research for detection in dynamical changes in electroencephalogram (EEG) signal with no consideration of the variability of this complexity measure. An application of the parametric bootstrap methodology is used to compare normal and pre-ictal EEG signals.

Tables

Table 1. Table 1: 90 % percent 90 90\% Confidence Intervals (Eq. 16 ) for every symbol length m 𝑚 m and power law parameter α 𝛼 \alpha . How many times the real value of ℋ ℋ {\mathcal{H}} - in fact we use the mean of ℋ ^ ​ ( T , m ) ^ ℋ 𝑇 𝑚 \hat{\mathcal{H}}(T,m) (see paragraph 4.1 ) that is the best possible estimator - is lower than the Lower Bound M ​ i ​ s ​ s ​ L ​ e ​ f ​ t 𝑀 𝑖 𝑠 𝑠 𝐿 𝑒 𝑓 𝑡 MissLeft or higher than the Upper Bound M ​ i ​ s ​ s ​ R ​ i ​ g ​ h ​ t 𝑀 𝑖 𝑠 𝑠 𝑅 𝑖 𝑔 ℎ 𝑡 MissRight . For white noise the confidence level is in fact higher than 90%, in fact is always accurate but for other values of α 𝛼 \alpha the overall confidence level is approximately between 90%. and 95%.

m	$\alpha$	$\hat{\mathscr{H}}(T,m)$	Miss left	MissRight	Mean Amplitude
3	-1	0.995831848	0	0.04	0.00222
4	-1	0.989083439	0.02	0.04	0.00420
5	-1	0.983495069	0.04	0.04	0.00500
6	-1	0.97547007	0.02	0	0.00555
3	0	0.99990292	0	0	0.00057
4	0	0.999679839	0	0	0.00080
5	0	0.998800463	0	0	0.00134
6	0	0.994503528	0	0	0.00235
3	1	0.991622896	0.02	0.02	0.00340
4	1	0.983385433	0.02	0.02	0.00493
5	1	0.97600538	0.02	0.04	0.00591
6	1	0.966355927	0	0.06	0.00657
3	2	0.943233315	0.08	0.06	0.00959
4	2	0.90634703	0.04	0.02	0.01273
5	2	0.878452628	0.02	0.02	0.01413
6	2	0.853327039	0.04	0.02	0.01482

Equations

$$\mathcal{H}[P]=S[P]/S_{max}=\left\{-\sum _{i=1}^N{P}_i\ln (P_i)\right\}/S_{max},$$

$$R(x_{t+n})=\sum _{k=0}^{m-1}\text{mathbbm}1(x_{t+k}\le x_{t+n})$$

$$X_t=(4,7,9,10,6,11,3)$$

$$P({\pi }_i)=\frac{{\sum }_{l=1}^{T-m+1}\text{mathbbm}1(X_m(l)\text{ has ordinal patter }{\pi }_i\text{ in }{S}_m)}{T-m+1},$$

$$\hat{\mathcal{H}}(m)=\left\{-\sum _{i=1}^{N}P({\pi }_i)\ln (P({\pi }_i))\right\}/S_{max},$$

$${\hat{P}}_T({\pi }_i)=\frac{n_i}{T-m+1}$$

$$P^{ij}=P(s_{t+1}={\pi }_j|s_t={\pi }_i)1\le i\le j\le m!$$

$$\hat{P}_{T}^{ij}=\begin{cases*}\frac{n_{ij}}{n_{i}}~{}~{}if~{}n_{i}\geq 0\\ 0~{}~{}~{}otherwise\end{cases*}$$

$$P({\pi }_j)=\sum _{i=1}^{m!}P({\pi }_i)P^{ij}$$

$${\hat{\mathcal{H}}}_T=\left\{-\sum _{i=1}^N{\hat{P}}_T({\pi }_i)\ln ({\hat{P}}_T({\pi }_i))\right\}/ln(m!)$$

$${\hat{\sigma }}_B({\hat{\mathcal{H}}}_T)=\hat{\sigma }({\hat{\mathcal{H}}}_T^{\ast })$$

$$\hat{\sigma }({\hat{\mathcal{H}}}_T^{\ast })=\sqrt{\frac{1}{B-1}\sum _{i=1}^B{\left({\hat{\mathcal{H}}}_T^{\ast }(i)-{\hat{\mathcal{H}}}_T^{\ast }(\bullet )\right)}^2}$$

$${\hat{\mathcal{H}}}_T^{\ast }(\bullet )=\frac{1}{B}\sum _{i=1}^B{\hat{\mathcal{H}}}_T^{\ast }(i)$$

$$\text{Bias}({\hat{\mathcal{H}}}_T^{\ast })={\hat{\mathcal{H}}}_T^{\ast }(\bullet )-{\hat{\mathcal{H}}}_T$$

$$\text{MSE}({\hat{\mathcal{H}}}_T^{\ast })=Var({\hat{\mathcal{H}}}_T^{\ast })+Bia{s}^2({\hat{\mathcal{H}}}_T^{\ast })$$

$$\left[\max (2.{\hat{\mathcal{H}}}_T-{\hat{\mathcal{H}}}_T^{\ast }(\bullet )+{\delta }_{\frac{\alpha }{2}}^{\ast },0),\min (2.{\hat{\mathcal{H}}}_T-{\hat{\mathcal{H}}}_T^{\ast }(\bullet )+{\delta }_{(1-\frac{\alpha }{2})}^{\ast },1)\right]$$

$${\mathcal{H}}_1\ne {\mathcal{H}}_2$$

$$\hat{\sigma }({\hat{\mathcal{H}}}_T)=\sqrt{\frac{1}{n-1}\sum _{i=1}^n{\left({\hat{\mathcal{H}}}_T(i)-{\hat{\mathcal{H}}}_T(\bullet \right)}^2}$$

$${\hat{\mathcal{H}}}_T(\bullet )=\frac{1}{n}\sum _{i=1}^n{\hat{\mathcal{H}}}_T(i)$$

Full text

††journal: Journal of LaTeX Templates

1 Introduction

In 2002, Bandt and Pompe introduced a measure of complexity for time series [4] named permutation entropy (PE). It is an information entropy [20] that takes account of the time evolution of the time series, in contrast with other prominent information entropies as the Shannon entropy [40]. It computation is fast, requires not too long time series [36] and it is robust against noise [33]. This measure has been widely used in non-linear dynamics [23, 15, 38, 26], and to a lesser extent in Stochastic Processes [37, 41, 46], among others. It has also had a great impact in in such different and important areas of applied science and engineering as varied as Mechanics Engineering [44, 34], Epilepsy [35, 31], anaesthesia [21], Cardiology [17, 32], Finance [28], Climate Change [12].

Since its publication and up to the end of 2016, this paper has been cited in 789 papers, according Scopus bibliographic database, and the evolution of the cites seems to indicate that it will be increasing within time. All these facts made an investigation of PE from the statistic point of view an important issue.

There has been little research, up to our knowledge, on the statistical properties of the quantities used in nonlinear dynamics to characterize time series. This lack of research may be due the lack of distributional theory of these quantitites, yielding resampling technics as the most powerfull tool to overcome this task. Perhaps one exception to this is the research on the distribution of the largest Lyapunov exponent and the correlation dimension[9]. We will summarize one of the most important discussion in this matter, according to our criterion and having in mind the computational scheme that we are proposing: In [18] a methodology to calculate the empirical distributions of Lyapunov exponents based on a a traditional bootstrapping technique is presented, providing a formal test of chaos under the null hypothesis. However, in [45] it is shown that the previously bootstrap approach seems to fail to provide reliable bounds for estimates the Lyapunvos exponents, and conclude that the traditional bootstrap cannot be applied for estimating multiplicative ergodic statistics. In [10], a moving blocks bootstrap procedure is used to detect a positive Lyapunov exponent in financial time series. However, the time series generated by moving block bootstrap present artifacts which are caused by joining randomly selected blocks, so the serial dependence is preserved within, but not between, the block.

Regarding time series symbolic dynamics, in [5] the probabilities generated using the Bandt and Pompe methodology are calculated analytically for Gaussian Processes for symbol length equal three, but they recognize that for larger length this is not possible, for that reason a computer based method is required to estimate the bias and variance in the PE estimation.

In this contribution we propose a different simulation method (i.e parametric bootstrap) for estimating the bias, variance and confidence intervals for the Permutation Entropy estimation, along with hypothesis testing, that consists in simulate bootstrap symbolic time series samples that are thought to be produced by a probabilistic model with a fixed transition probability extracted from the original time series.

In order to show some results from our method we simulate a well known family of time series: the $1/f^{\alpha}$ noise. We compute bias, variance and confidence intervals for the Permutation Entropy of these time series according time series length and several parameters. In addition, an application of the parametric bootstrap methodology for hypothesis testing is used to compare normal and pre-ictal EEG signals.

The paper reads as follows: Section 2 shows a brief review of PE in order to present the estimator to be evaluated using the bootstrap approach, Section 3 presents and explains the proposed parametric bootstrap, firstly a brief review of the bootstrap scheme is done as introduction to our method, then in Subsection 3.1 the core of the bootstrap approach is presented, i.e. the probability transitions computation is explaiedn and finally in Subsection 3.2 the algortihm to parametric bootstrap PE is explained. Section 4 presents the dynamical systems simulated, Section 5 introduces the experimental data used in the application and Section 6 is devoted to the results and conclusions of this contribution.

2 Permutation Entropy

In this Section we briefly review the PE to make the article self contained and accessible for a wider audience.

Let $\{X_{t}\}_{t\in T}$ be a realization of a data generator process in form of a real valued time series of length $T\in\mathbb{N}$. A measure of uncertaintly about $\{X_{t}\}_{t\in T}$ is the normalized Shannon entropy [40] ($0\leq{\mathcal{H}}\leq 1$), which is defined as:

[TABLE]

where $P_{i}$ is a probability to be extracted from the time series, $N$ is the cardinality of the $P_{i}$ set $\{p_{i}\}_{1}^{N}$, the denominator $S_{max}=S[P_{e}]=\ln N$ is obtained by a uniform probability distribution $P_{e}=\{P_{i}=1/N,~{}\forall i=1,\cdots,N\}$.

Bandt and Pompe proposed a symbolization technique to estimate $P_{i}$ and compute PE, $\hat{\mathcal{H}}(m,\tau)$. First we recall that PE has two tuning parameters, i.e. $m$ the symbol length and $\tau$ the time delay. Within this paper, we set $\tau=1$ with no loss of generality and it will omitted, so we will use $\mathcal{H}=\mathcal{H}(m)$ for sake of simplicity. Let $X_{m}(t)=(x_{t},x_{t+1},\dots,x_{t+m-1})$ with $0\leq t\leq T-m+1$ be a non-disjoint partition containing the vectors of real values of length $m$ of the time series $\{X_{t}\}_{t\in T}$. Let $S_{m\geq 3}$ the symmetric group of order $m!$ form by all possible permutation of order $m$, $\pi_{i}=(i_{1},i_{2},\dots,i_{m})\in~{}S_{m}$ ($i_{j}\neq i_{k}\forall j\neq k$ so every element in $\pi_{i}$ is unique). We will call an element $\pi_{i}$ in $S_{m}$ a symbol or a motiv as well. Then $X_{m}(t)$ can be mapped to a symbol $\pi_{i}$ in $S_{m}$ for a given but otherwise arbitrary $t$. The $m$ number of real values $X_{m}(t)=(x_{t},x_{t+1},\dots,x_{t+m-1})$ are mapped onto their rank. The rank function is defined as:

[TABLE]

where $\mathbbm{1}$ is the indicator function (i.e $\mathbbm{1}(Z)=1$ if $Z$ is true and [math] otherwise) , $x_{t+n}\in X_{m}(t)$ with $0\leq n\leq m-1$ and $1\leq R(x_{t+n})\leq m$. So the rank $R(min(x_{t+k}))=1$ and $R(max(x_{t+k}))=m$. The complete alphabet is all the possible permutation of the ranks. Hence, any vector $X_{m}(t)$ is uniquely mapped onto $\pi_{i}=\left(R(x_{t}),R(x_{t+1}),\dots,R(x_{t+m-1})\right)\in S_{m}$. With this Rank Permutation Mapping one simply maps each value $x_{i}$ in $X_{m}(t)$ placing its rank $R(x_{i})\in\{1,2,\dots,m\}$ in chronological order to form $\pi_{i}$ in $S_{m}$. In Figure 1 an illustrative drawing of this mapping for all alternatives in $m=3$ is presented. It can be seen that the indexes of the vertical axis are fixed, ordered by amplitude (i.e ranks), and they are mapped onto the time axis. The resultant symbol can be obtained reading the labels in the horizontal axis from left to right (in chronological order). This method is used by [3, 36, 6] among others. For example, let us take the series with seven values ($T=7$) [4] (see Fig. 2, top), and motiv length $m=3$:

[TABLE]

$X_{3}(1)=(4,7,9)$ and $X_{3}(2)=(7,9,10)$ represents the permutation $\pi=123$ since $R(x_{1})=1$ ,$R(x_{2})=2$, $R(x_{3})=3$. $X_{3}(3)=(9,10,6)$ and $X_{3}(4)=(6,11,3)$ correspond to the permutation $\pi=231$ since $R(x_{1})=2$ ,$R(x_{2})=3$, $R(x_{3})=1$ (see Fig. 2, middle). Using the rank permutation Mapping we compute $P(\pi_{i})$ (see Fig. 2, bottom) ,

[TABLE]

where $\mathbbm{1}$ is the indicator function and $i=1,\dots,m!$. Using these probabilities, $\hat{\mathcal{H}}(m)$ can be computed as,

[TABLE]

where $N=m!$ is the order of the symmetric group $S_{m}$ and $S_{max}=log(N)$.

3 The Bootstrap approach

The bootstrap is a computer based method for assigning measures of accuracy to the desired statistical variable estimates. If $\mathcal{H}$ is an unknown characteristic of a model $\Psi$, an estimator $\hat{\mathcal{H}}$ can be derived from the sample generated by $\Psi$ in a single experiment. A way to obtain the distribution of $\hat{\mathcal{H}}$ is to repeat the experiment a large number of times and approximate the distribution of $\hat{\mathcal{H}}$ by the so obtained empirical distribution. In most practical situations this method is impossible because the experiment is not reproducible, or is unenforceable for cost reasons. The spirit of the bootstrap methodology is to estimate the sampling distribution of a statistic (i.e quantifier or a parameter estimator) from the data at hand by analogy to the ’thought experiment” that motivates the sampling distribution.

Suppose that an unknown probability model $\Psi$ gives an observed data set $\textbf{X}=\{x_{t}\}_{t=1}^{n}$ by random sampling and let $\hat{\theta}_{\Psi}(\textbf{X},T)$ be the statistic of interest that estimates our true value $\theta=f(\Psi)$. Then with the observed data set X we produce an estimate $\hat{\Psi}$. The trick is now to repeat the random sampling but with the estimate $\hat{\Psi}$ giving bootstrap samples $\textbf{X}^{*}=\{x^{*}_{t}\}_{t=1}^{n}$ and for each bootstrap sample we calculate $\hat{\theta}^{*}_{\hat{\Psi}}(\textbf{X}^{*},n)$. Now we repeat the bootstrap sampling $B$ times and the distribution of $\hat{\theta}^{*}_{\hat{\Psi}}(\textbf{X}^{*},n)$ is the bootstrap estimator of the distribution of $\hat{\theta}_{\Psi}(\textbf{X},n)$. With this estimated distribution we can also obtain the variance, the bias and the confidence intervals of our estimator.

Formally, the bootstrap methodology is based in the plug-in principle. The parameter of interest can be written as a function of the probability model, $\theta=f(\Psi)$. As the probability model is unknown, the plug-in estimate of our parameter is defined to be $\hat{\theta}=f(\hat{\Psi})$. So the bootstrap propose that we resample from this estimated probability model $\hat{\Psi}$ that is chosen to be close to $\Psi$ in some sense.

If we have some information about $\Psi$ besides the data, the chosen $\hat{\Psi}$ must contain this information. Suppose that we know that the data $\textbf{X}=\{x_{t}\}_{t=1}^{n}$ comes from a certain process ruled by a probabilistic model $\Psi$ that depends on a finite number of parameters $\mathbf{\Xi}=\{\xi\}_{i=1}^{k}$, so $\Psi~{}=~{}f\left(\mathbf{\Xi}\right)$. This parameters can be estimated in the traditional statistical parametric approach as Maximum Likelihood getting $\mathbf{\hat{\Xi}}=\{\hat{\xi}\}_{i=1}^{k}$, and equivalently $\hat{\Psi}~{}=~{}f\left(\mathbf{\hat{\xi}}\right)$.

Now the bootstrap samples $\textbf{X}^{*}=\{x_{1}^{*},x_{2}^{*}\dots x_{n}^{*}\}$ comes from a process ruled by a probabilistic model $\hat{\Psi}$ (Fig. 3). These bootstrap samples emulates in every sense the original samples, including the correlation between the values.

3.1 The transition probabilities of a symbol sequence

As stated in Section 2 using the methodology proposed by Bandt & Pompe, the dynamics of a process $\{X_{t}\}_{t\in T}$ with $X_{t}\in\mathbb{R}$ is represented by a $m!-th$ finite state random process $\{S_{t}\}_{t\in(T-m+1)}$ with $S_{t}\in S_{m}=\{\pi_{1},\pi_{2}\dots\pi_{m!}\}$ for all posible $m\geq 2$. This realization of the symbolic sequence is thought to be produced by a probabilistic model with a fixed transition probability $P^{ij}$ (i.e. the probability of moving from a symbol $\pi_{i}$ to a symbol $\pi_{j}$ for all $1\leq i\leq m!$ and $1\leq j\leq m!$) denoted $\Psi(P^{ij})$.

We estimate the model $\hat{\Psi}=\Psi(\hat{P^{ij}})$, see Fig. 3 to bootstrap the Permutation Entropy

According to [4], the relative frequency

[TABLE]

is an estimate as good as possible for a finite series of values of $P(\pi_{i})$, with $n_{i}$ the number of times the state $\pi_{i}$ is observed up to time $T-m+1$. The sub-index $T$ in $\hat{P}_{T}(\pi_{i})$ reinforces the notion of the dependence of the estimator on the length of the series $T$.

With the same spirit we define the transition probabilities of the symbol sequence as:

[TABLE]

And the estimator of $P^{ij}$

[TABLE]

where $n_{ij}$ is the number of transitions observed form $\pi_{i}$ to $\pi_{j}$ up to time $T-m+1$ and . Note $n_{i}~{}=~{}\hat{P}(\pi_{i}).(T-m+1)$

Then by the law of the total probability:

[TABLE]

so if we call $\mathbf{P}(\pi)$ to the (m!)-dimensional vector containing $P(\pi_{i})$ in each coordinate (i.e $\mathbb{P}(\pi)~{}=(P(\pi_{1}),P(\pi_{2}),\dots,P(\pi_{m!})$), then $\mathbf{P}(\pi)$ is determined by $P^{ij}$, leading to the conclusion that the estimator $\hat{\mathbf{P}}_{T}(\pi)$ is determined by the estimation of $\hat{P}_{T}^{ij}$.

3.2 Bootstrapping the Permutation Entropy

The Permutation Entropy is defined in eq. 1, so because of the plug-in principle, our natural estimator is:

[TABLE]

In section 3.1 we showed that the Permutation Entropy was completely defined by the transition probabilities $P^{ij}$ so we can think of them as parameters of a probabilistic model $\Psi$.

Following the scheme in Fig. 3 we have:

$\Psi\left(P^{ij}\right)\longrightarrow\mathbf{S}=(s_{1},s_{2}\dots,s_{T-m+1})\longrightarrow\hat{\mathcal{H}}_{T}$

Our probabilistic model with unknown transition probabilities $P^{ij}$ gives the observed symbol sequence $\mathbf{S}$, and with that sequence the estimation of the Permutation Entropy is obtained.

In the ’bootstrap world”:

$\hat{\Psi}~{}=~{}\Psi\left(\hat{P}_{T}^{ij}\right)\longrightarrow\mathbf{S^{*}}=(s_{1}^{*},s_{2}^{*}\dots,s^{*}_{T-m+1})\longrightarrow\hat{\mathcal{H}}^{*}_{T}$

$\hat{\Psi}$ generates $\mathbf{S^{*}}$ by a simulation, giving the bootstrap replication $\hat{\mathcal{H}}^{*}_{T}$. We can repeat the simulation to get as many bootstrap replications as affordable.

Computing $B$ bootstrap replication of the permutation entropy from a time series $\{X_{t}\}_{t\in T}$ is simple: given a time series of lenght $T$, choose a world length $m$ and a time delay $\tau$ to do the mapping from $\{X_{t}\}_{t\in T}$ to $\{S_{t}\}_{t\in(T-m+1)}$ as stated in section 2. With this sequence: compute $\hat{P}_{T}(\pi_{i})$, (eq. 6), $\hat{P}_{T}^{ij}$ (eq. 8) and calculate $\hat{\mathcal{H}}_{T}$ (eq. 10). Then choose at random with probability $\hat{P}_{T}(\pi)$ an inicial state $s^{*}_{1}(b)=\pi_{k}$ and choose at random with probability $\hat{P}_{T}^{kj}$ (note that $k$ is fixed with the value of the previous state) the next simulated state $s_{2}^{*}(b)$. Repeat this last step $T-m+1$ times to obtain the simulation $\mathbf{S^{*}(1)}=(s_{1}^{*}(1),s_{2}^{*}(1)\dots,s^{*}_{T-m+1}(1))$. With this bootstrap replication of symbol sequence estimate $\hat{\mathcal{H}}^{*}_{T}(b)$ (equation 6).

For a more detailed reference see Algorithm 1 in Appendix A.

Repeat the simulation of the sequence $B$ times to obtain $\hat{H}^{*}_{T}(b)~{}~{}b=1\dots B$. With the set $\hat{\mathcal{H}}^{*}_{T}(b)~{}~{}b=1\dots B$ we have the bootstrap replications needed to estimate the standard deviation, the confidence intervals of $\hat{\mathcal{H}}_{T}$, or the test presented in the following section.

So we obtained $B$ bootstrap replications of $\hat{\mathcal{H}}^{*}_{T}$:

$\hat{\mathcal{H}}^{*}_{T}(1),\hat{\mathcal{H}}^{*}_{T}(2)\dots\hat{\mathcal{H}}^{*}_{T}(B)$

The Bootstrap Standard Deviation of $\hat{\mathcal{H}}^{*}_{T}$ is our estimation of the Standard Deviation of $\hat{\mathcal{H}}_{T}$:

[TABLE]

and is defined as

[TABLE]

where

[TABLE]

We define the bootstrap bias of $\hat{\mathcal{H}}^{*}_{T}$ as:

[TABLE]

Finally, the Mean Square Error (MSE) of an estimator:

[TABLE]

3.2.1 Confidence Intervals

The $1-\alpha$ Confidence Interval of $\mathcal{H}$ is defined by the percentiles of the bootstrap $\delta$. For each bootstrap replicate $\hat{\mathcal{H}}^{*}_{T}(b)$ we compute the difference - $\delta^{*}(b)$ - between that replication and the mean of all bootstrap replicates. Then we choose the ($\frac{\alpha}{2}$) and the ($1-\frac{\alpha}{2}$) percentiles of the $\delta^{*}$’s distribution and add them to the original estimate, $\hat{\mathcal{H}}^{*}_{T}$, correcting for the bias, and the resulting $(1-\alpha)100\%$ confidence interval is:

[TABLE]

For a more detailed reference see Algorithm 2 in Appendix A.

3.2.2 Hypothesis testing

With this same spirit, a confidence interval for the difference between the permutation entropy of two different time series can be made. In inferential statistics exists a direct relationship between confidence intervals and hypothesis testing. A two-sided $(1-\alpha)$ confidence interval in the difference between two measures can be used to determine if those two measures are significantly different by only checking if the $zero$ belongs to this particular interval.

$H_{0}:\Delta=\mathcal{H}_{1}-\mathcal{H}_{2}=0$

If $0~{}\notin~{}(1-\alpha)100\%~{}CI~{}(\Delta)$

then reject $H_{0}$ and

[TABLE]

The procedure to perform this test is shown in 3 in Appendix A.

4 Numerical simulation

In order to show our proposed bootstrap in a very general time series, we simulate a well known dynamical system: the $1/f^{\alpha}$. All the series are simulated with different time spam, $T$ in order to evaluate the statistical properties of $\hat{\mathcal{H}}_{T}$ according to Ec. 10. As stated before, a way to obtain the distribution of $\hat{\mathcal{H}}$ is to repeat the experiment a large number of times and approximate the distribution of $\hat{\mathcal{H}}$ by the so obtained empirical distribution. While for real world experiments this can be inapplicable, for simulated time series this can easily done by Montecarlo Simulation. Once the $n$ replications of $\hat{\mathcal{H}_{T}}~{}=~{}\{\hat{\mathcal{H}_{T}}(1),\dots,\hat{\mathcal{H}_{T}}(n)\}$ is obtained the standard deviation is estimated by:

[TABLE]

where

[TABLE]

4.1 Experimental design

Stochastic dynamical systems:

$1/f^{\alpha}$ noises refers to a signal with spectral density $S(f)$ with the form $S(f)~{}=~{}k\frac{1}{f^{\alpha}}$ where $k$ is a constant, $\alpha$ is the signal-dependent parameter and $f$ is frecuency [22]. It is a stochastic model which seems to be ubiquitous in nature [22] and the references therein. We simulate $1/f^{\alpha}$ noises with $\alpha=\{-1,0,1,2\}$. See Fig. 4 for an example of these noises. A white noise process ($\alpha=0$) would generate a curve with constant power in the spectrum. The case of $\alpha=1$ or pink noise is the canonical case and of most interest as many of the values of $\alpha$ found in nature are very near to 1.0 [11, 24, 29, 43, 16]. A random walk noise (Brownian Motion or red noise, $\alpha=2$) would show a $(1/f^{2})$ distribution in $S(f)$. In order to simulate this stochastic process, the algorithm propose in [42] is used.

For each $\alpha=\{-1,0,1,2\}$ 1000 replications were simulated for each $T=\{60,100,120,400,$ $600,2000,3600,5000,10000,20000,50000\}$.

$\hat{\mathcal{H}}_{i}(T,m)$ for $\{i=1\dots 1000\}$ along with $\hat{\sigma}(\hat{\mathcal{H}}_{T})$ are obtained for $m=\{3,4,5,6\}$.

For each $\alpha=\{-1,0,1,2\}$ a single replication was simulated for each $T=\{60,100,120,400,600,2000,3600,5000,$ $10000,20000,50000\}$. In each case for this replication we implemented the algorithm 1 to get 1000 bootstrap replicates. $\hat{\mathcal{H}}_{i}^{*}(T,m)$ for $\{b=1\dots 1000\}$ and $\hat{\sigma}_{B}(\hat{\mathcal{H}}_{T})$ are obtained for $m=\{3,4,5,6\}$. For these bootstrap distribution we analyze Bias, Standard Deviation and MSE.

As for each set of 1000 bootstrap replicates we obtain a single confidence interval, we repeated this step 50 times to obtain Table 1 that indicates the estimated confidence level of this method along with the mean amplitude of the interval.

5 Application: EEG data

In order to illustrate the proposed confidence intervals in real contexts we present how it can describe the variability in the Permutation Entropy within one observation of Electroencephalogram (EEG) Data. More precisely, as a first practical application, we analyze, via PME, four different sets of EEGs for healthy and epileptic patients that were previously analyzed by [1]

(available at http://www.meb.unibonn.de/epileptologie/science/physik/eegdata.html). The data consist of 100 data segments (from which we choose 10 at random), whose length is 4097 data points with a sampling frequency of 173.61Hz, of brain activity for different groups and recording regions: surface EEG recordings from five healthy volunteers in an awake state with eyes open (Set A) and closed (Set B), intracranial EEG recordings from five epilepsy patients during the seizure free interval from outside (Set C) and from within (Set D) the seizure generating area. Details about the recording technique of these EEG data can be found in the original paper.

6 Results and discussion

We intend to show in our simulated experiment that the bootstrap distribution of the PE estimator is close in every meaningful sense to the distribution obtaining by the repetition of the original experiment (empirical distribution) in order to obtain this estimator distribution when the exact replication of the experiment can not be done. A comparison between the standard deviations of both bootstrap replicates and the empirical distribution ($\hat{\sigma}_{B}(\hat{\mathcal{H}}_{T})$ and $\hat{\sigma}(\hat{\mathcal{H}}_{T})$ respectively) for the stochastic processes is presented in Fig. 5. There are some discrepancies for low values of $T$, but form a certain value $T_{0}$ in all the cases of $m$ the standard deviation coincides. In Fig. 6 it can be seen that for every $m$ and $\alpha$ the bias of the bootstrap estimate tends to zero as $T$ increases. So, this bootstrap estimator is an asymptotically unbiased estimator. With this and with the fact that $\sigma$ also goes to zero as $T$ increases, the bootstrap estimator seems to be Mean Square Consistent. Even more, for large values of $T$, the bootstrap estimation is as efficient as the estimation produced by the repetition of the experiment.

In Fig. 7 and for an arbitrary value of $\alpha=1$ an for the largest length of the simulated time series $T=50000$, an histogram of the bootstrap estimator along with an histogram of the simulated estimator are presented in different scale for every $m$. The similar shape between the histograms can be appreciated, the difference in the location is due that the bootstrap samples depends on only one of the estimations of the PE (that are random) but this does not affect the posterior inferential conclusions.

For a more thorough exploration of the bootstrap estimator we have calculated fifty $90\%$ Confidence Intervals and for every $m$ and every $\alpha$ we computed how many times the real value of ${\mathcal{H}}$ - in fact we use the mean of $\hat{\mathcal{H}}(T,m)$ (see paragraph 4.1) that is the best possible estimator - is outside the bounds of the confidence interval. Results are shown in Table 1. For white noise the confidence level is in fact higher than 90%, in fact is always accurate but for other values of $\alpha$ the overall confidence level is approximately between 90%. and 95%.

In many practical situations, there is a wish to compare the dynamics of two processes via the Permutation Entropy of their time series. The question is: $\mathcal{H}_{1}=\mathcal{H}_{2}$? This can not be answered with punctual estimators ($\hat{\mathcal{H}}_{1},\hat{\mathcal{H}}_{2}$) because these are continuous random variables and with probability 1 (i.e. always) they are going to be different. The real question is if that difference is statistically significant or not, and that only can be answered if exists a measure of variability of that continuous random variable, $\hat{\Delta}=\hat{\mathcal{H}}_{1}-\hat{\mathcal{H}}_{2}$. There has not been, up to our best knowledge, this kind of variability measure that we are proposing now. An example of this is the Permutation Entropy of EEG signals.

The problem of interest is comparing the PE of 4 different sets of EEG signals: EEG signals of patients in an awake state with eyes open (Set A) and closed (Set B), intracranial EEG recordings from epilepsy patients during the seizure free interval from outside (Set C) and from within (Set D) the seizure generating area.

If many EEG signals for each type of patients can be recorded a classical inference for the mean can be performed if the normality assumptions are complied or, if normality fails (that is to be expected in this case), a non parametric test can be made. But is the mean PE representative of the population of each type of patient?

A different problem is to analyze the variability of a single EEG signal, this can not be done with conventional methods and up to this contribution there has not been an answer to this problem. The same problematic applies when the issue is to compare between two EEG signals.

We solve this problem by constructing confidence intervals and hypothesis testing with our proposed method. In Fig 8 90% Confidence Intervals for the 10 EEG signals of brain activity for different groups and recording regions are performed. It should be pointed out that the overlapping between intervals does not necessarily means that there is no significant differences between the two Permutation Entropies. To reach that conclusion, an hypothesis test for the difference must be made.

In Fig. 9 we perform a test for difference in the Permutation Entropy between the 10 EEG signals of healthy volunteers in an awake state with eyes open ($SetA$) and the 10 EEG signals of healthy volunteers in an awake state with eyes closed ($SetB$). Each EEG signal of $SetA$ was compared with each signal of $SetB$ with a 10% significance level, and the conclusion is that the differences seems to be at random, indicating that is no real difference between these two types of EEG signals. In Fig. 10 the same analysis is extended to all the different types of patients. While the differences between $SetA$ and $SetB$ seems to be at random, all EEG signals of those Sets are different in every test to the EEG signals of $SetC$ and $SetD$. Instead, between $SetC$ and $SetD$ again the differences are distributed between significant and not significant.

In summary, we present a computer based methodology to obtain an accuracy measure for the estimation of the permutation entropy ($\hat{\mathcal{H}}$). So far we found in the literature that only descriptive statistics are used to characterize this quantifier and if the objective is to extrapolate on and reach conclusions that extend beyond the raw data itself there were no statistical inference method at hand. Even a simple comparison between to random variables (as $\hat{\mathcal{H}}$) can not be made with some confidence without a measure of variability of that variable. Our method paves the way to perform any inferential statistic involving the Permutation Entropy or even any entropy that uses the Probability Function Distribution proposed by Bandt and Pompe.