Five snippets of signal analysis

Signal Analysis currently features tests and procedures addressed at verifying certain hypotheses, such as linearity, stationarity, changes of regime, and the kind of information utilized in back/forecasting time series, possibly by calibration. In such context, methodological updates are required, especially when it comes to dealing with series characterized by high irregularities or when unobservable are involved as in the case of Kalman filtering.1,2 Many Regime-Switching techniques, for instance, are affected by strict requirements such as stationarity and single-point shift3,4 which require amelioration or updating. Similarly, time-varying parameters utilized for back/ forecasting purposes are better computed in a Bayesian context, in order to exploit prior pieces of evidence based on their known distribution. Finally, current stationarity tests5,6 in a nonlinear setting with dramatic changes in trend and subsample variances may require substantial corrections. Section 1 is devoted to defining random signals, whereas Section 2 tackles the issue of Regime-Switch testing. Section 3 produces a novel technique in the field of signal smoothing. The following two sections illustrate the steps for implementing Bayesian Monte Carlo simulation and calibration, whereas the final section is addressed at supplying corrected stationarity tests. Section 7 concludes.

Introduction
Signal Analysis currently features tests and procedures addressed at verifying certain hypotheses, such as linearity, stationarity, changes of regime, and the kind of information utilized in back/forecasting time series, possibly by calibration. In such context, methodological updates are required, especially when it comes to dealing with series characterized by high irregularities or when unobservable are involved as in the case of Kalman filtering. 1,2 Many Regime-Switching techniques, for instance, are affected by strict requirements such as stationarity and single-point shift 3,4 which require amelioration or updating. Similarly, time-varying parameters utilized for back/ forecasting purposes are better computed in a Bayesian context, in order to exploit prior pieces of evidence based on their known distribution. Finally, current stationarity tests 5,6 in a nonlinear setting with dramatic changes in trend and subsample variances may require substantial corrections. Section 1 is devoted to defining random signals, whereas Section 2 tackles the issue of Regime-Switch testing. Section 3 produces a novel technique in the field of signal smoothing. The following two sections illustrate the steps for implementing Bayesian Monte Carlo simulation and calibration, whereas the final section is addressed at supplying corrected stationarity tests. Section 7 concludes.

Definition of signal
All of the aforementioned procedures share the assumption of a signal represented by a real-valued random sequence which, for the discrete time notation [ ] With unknown distribution, mean 1 1 T X T X t t − ∑ = = and finite variance 2 x σ < ∞ . More specifically, Equation (1) may be decomposed Cramér 7 in the following fashion * X X S t t t t ε = + + Where ( )

Sequential and cluster-based regime switching
Regime switching is intended as a single or multiple change of pattern of the signal in terms of mean, variance, and autocorrelation structure. 9,10 It is often associated to rare events such as disasters with permanent effects like earthquake aftermaths or financial crises. 11 Changes associated to regime switching are generally abrupt and of sizable magnitude, and may be tested for their statistical significance by means of some preselect critical-value standard. 12 Regime switching is closely associated to structural-break testing. 13 However, it is more pervasive by attempting to establish not only the date and size of breaks, but also the moment effects over the signal of interest. Moreover, regime switching is deeply involved with threshold behavior, which analyzes the event of a variable taking an extreme value at some time of its evolution. In practice, the parameterization of the signal after each break is tested for in a context of a Markov switching model. Many tests exist for the detection of regime switches, 2 whereby the different states of the signal may be identified within a given degree of likelihood. Among these tests, two stand out for their usefulness and simplicity: Sequential F-test (two-sample Kolmogorov-Smirnov test) for linear (nonlinear) signals. 3,4 Sample cluster-based Quasi-Likelihood Ratio (QLR) test of the τ∈ Τ Τ< . The subsample length is established arbitrarily, usually on a peak-to-peak basis. The F-test version for establishing the number of regime switches is given by the following statistic

Var
represents the variance of the braced argument. From Equation (3), the M-sized string of F statistics and associated p-values may be obtained. A regime switch, namely a subsample variance outlier with respect to the sample variance, is said to significantly exist if its m.th p-value falls short of some preselect critical value, usually 5%.
In case of a nonlinear signal, the test statistic corresponding to F m is represented by the two-sample Ansari-Bradley test for equal sample tied-rank scores, also defined as quasi-variances. 15 The Ansari-Bradley test requires several steps departing from the padded signal vector { ; }; = + ; = + ; , üüü m t θ τ θ τ = ∈Θ Θ Τ ∀ ∈ , and from the associated vector R θ of the locations of the values in Z θ ranked in increasing order. Next, the locations of vector R θ are parted into two group vectors depending on their origin, namely, their belonging to the subsample or to the entire sample. Formally, the concept may be expressed as Subsequently, the tied-rank vector U θ from each end of Z θ is computed, so that the smallest and largest values of U θ get rank 1, the next smallest and largest get rank 2, etc. Finally, the Ansari-Bradley test statistic 15 is shown to be computed as follows Which is the sum of the tied ranks associated to the m.th subsample? From Equation (4), the test statistic * W m may be obtained by subtraction of the mean tied-rank score sum together with the corresponding p-value. The statistic is asymptotic normal under the null hypothesis that two stochastic processes come from the same distribution, against the alternative that they come from distributions sharing the same median but different variances. In the present context, it may also be loosely viewed as a test for heteroskedasticity. 16 The QLR test, in its simplest form, checks the null hypothesis of one versus two regime switches, and is built up from forming two non-sequential clusters each characterized, possibly, by different distributions. Thereafter, the QLR is computed upon the two samples, with the proviso that these produce a Gaussian mixture where the standard Likelihood-Ratio testing cannot be implemented. For this purpose, the relevant literature provides appropriate critical values. 12 Because of its stricture in admitting a maximum of two regime switches and of its computational complexity, the QLR test receives here no applied treatment.

Spectral representation theorem for signal decomposition and optimal smoothing
A large amount of signal smoothers (i.e. smoothed trends) is nowadays available for applied work in many different scientific fields. Of these, the most notable are the Savitzky-Golay, the Daubechies, the Hodrick-Prescott filters and the Hilbert Huang Transform. [17][18][19][20][21][22][23] All of these methods, however, share serious problems when it comes to data resolution. In fact, they all require prior input of specific parameters, namely, the polynomial degree for the first two methods, the window width for the first and the decomposition level for the second method. In addition, it requires a preselected smoothing coefficient for the third method and, finally, a 'sifting' grid search for the last.
For details of such inputs, the reader is redirected to the respective reference listed in this paper. Suffice here to warn that arbitrary, and thus casual, choice of any of these parameters is conducive to the risk of over/under fitting of the slow mode * X t with respect to the original signal X t of Equation (1). The immediate consequence is inefficient estimation in terms of suboptimal denoising, namely and loosely stated, a bang-bang signal extraction yielding either a too flat smoother that is statistically inconsistent with the original peak/ trough pattern or that too closely replicates the original signal itself. We propose here a fast, automatic, and efficient smoothing procedure that applies the Spectral Representation Theorem to a signal 1 X of whichever nature: stationary or not, linear or nonlinear. The signal may be approximated by means of the following periodic function where ˆ, X t k is a fitted sequence of the original signal, µ is the is a functional form, ( ) . X f t ∆ = expresses the existence of linearity and nonlinearity of the process as evidenced by the value of the exponent attached to the lagged signal. The fitted signal of Equation (5), which is the smoother estimable by Ordinary Least Squares (OLS), may be defined as the Spectral Representation Transform of the original signal with time-varying amplitude (Hamilton, 1994 . The optimal smoother is therefore defined as For the optimal lag * K of Equation (5) and the performance index closely akin to the Percent Root Mean Square Error (PRMSE), defined as Where . is the Euclidean norm of its argument.

Bayesian monte carlo simulation with bootstrapping for back/forecasting
Bayesian Monte Carlo simulation with bootstrapping is one of the two signal back/forecasting procedures herein proposed. The other procedure is shown in Section 6. Whereas the latter is a full-blown Kalman filter with time-varying parameters extracted from Gibbs sampling, 2 this procedure is based on bootstrapping known dates and values of the signal to produce a feasible MCMC which represents the back/forecasted signal along select P<∞ periods. This procedure is thus a game of chance which involves all of the elements of Equation (2) such that we may proceed to back/forecasting a signal by means of two different techniques: I. Back/forecast the smoother of Equation (4) by P periods and thereafter append stationary shocks drawn from the distribution of the original signal, 24 II. Perform random reconstruction of the cyclical pattern of the original signal by applying the principles of peak/trough selection and bootstrapping these cyclical extrema without replacement. 25 While the first technique works well only under the assumption of fully or fairly regular cycles in amplitude and length, and is likely to be affected by substantial uncertainties, 26-28 the second technique preserves the stochastics of the original signal and is by consequence more trustworthy from the distributional viewpoint. There follows that Bayesian prior values drawn from the known distribution of the original signal can be utilized along the entire signal sample so as to forestall predictive uncertainties. Implementation of this preferred technique requires three consecutive steps: Find the peaks and troughs of the signal X t throughout the entire sample. Scrap all occurrences separated by a cyclical distance far exceeding or falling too short of the average cycle.
The number of trough-to-trough cycles remained may be denoted as1 C T ≤  . Subsequently form a matrix ( ) ,3 F C of triplets of dates in the following order: lower-date trough, peak date, and upper-date trough.
The latter is the lower-date trough of the subsequent cycle, and so on. Each triplet encapsulates a specific cycle of given length. Bootstrap N times the rows of matrix F so as to obtain a matrix ( ) , G P D of the reshuffled values of the signal X t . For each row of this matrix such values might be further simulated up to D draws to produce Bayesian priors from the distribution of X t and generate a new matrix G with the posteriors. From this matrix find the series which mostly approximates the original signal X t in terms of mean and/or variance. This is the optimal prediction signal ** X p , for [ ] 1, p P ∈ , commonly defined as the "mode" since it is the closest to the actual data. This procedure is quick and easy to implement as opposed to the lengthier procedure featured in the next section, which needs a larger, often cumbersome informational set provided by the proxy instrument (s) utilized.

Bayesian calibration for back/forecasting
Calibration is a technique for back/forecasting the signal X t by means of a matrix of proxy instruments ( ) where T>T τ∈ and 1 n≥ respectively are the sample length and size. It assumes a statistically significant relationship among the two datasets for the overlapping sample length T such that back/forecasting P periods ahead (behind) by means of Z p is expected to produce a mode ** X p . The relationship is customarily expressed in terms of a fixed-time parameter set ( ) : ,1 0 B n usually estimated by means of OLS. 29 Bayesian calibration in the present context is an extension of the above model as it contemplates time-variable parameters connecting the two datasets and letting the prediction process be estimated by a dynamic statespace (SS) model whose upshot is a MCMC process characterized by ergodicity. 2 This model is represented by two equations whose workings in a Bayesian framework are better detailed elsewhere 30 Where, for p given as above (Sect. 1.3), and B X p p respectively are the state(s) and the observable(s), whereas and e v p p are orthogonal I.I.D. zero-mean disturbances. More precisely,

( )
: , B P n p is the time-varying parameter set whose first input is given by the Bayesian prior 0 B conditional upon a normal distribution. The variable C p is a slope parameter set connecting the contemporary state(s) and observable(s) and has the same size as B p . , n n except the last which is a scalar. The first component of p Ω is the Riccati matrix of the state covariance, which is distributed as an Inverse Wishart, whereas the other two components respectively are the state and the observable error covariances, distributed as Inverse Gamma. Matrix p Ω may be defined as the prior information set for predicting 1 B p+ and its estimation errors e p . Once this value is computed conditional on the covariances and on their distributional properties, the ensuing posterior 1 p Ω + is utilized for computing 2 B p+ in the same fashion, and so on along the process. This prior-posterior-prior sequence of the information set is at the heart of the conditional adaptive Kalman Filter for which the conditional parameters of Equation (9) at any stage are optimally estimated and ends as soon as P is hit. 1 Needless to say, optimality refers in the present context to the valuable information set p Ω delivered stage-wise onto the parameters. More explicitly, let the Kalman Gain be ( ) Such that the Riccati matrix be represented overtime as the following nonlinear and recursive process Needless to say, all of the elements in Equations (10) to (12) enter the real-time information matrix p Ω . Optionally, in order to gain further information by letting and 1 B B p p+ be conditioned directly by the observable X p , the covariance matrix Q p in Equation (11) may be replaced by the variance of the observable itself.

Proposed method for correction of stationarity test statistics
In this last Section a novel method is proposed for correcting the conventional stationarity test statistic in a setting of non-flat linear trends, and non-constant smoothers and cycles. Among such tests, the most notorious applicable to linear signals is the t-type Augmented Dickey-Fuller (ADF) test 5,6 which embeds an option for the intercept and/or for linear trend. Similarly, for nonlinear signals the Kapeitanos-Shin-Snell (KSS) t-type test is commonly applied with similar options. 31,32 The critical asymptotic values with preselected p-values of the test statistics are tabulated in the respective references. The null hypothesis implied by both test statistics involves first differencing of the signal, which in Econometrics is dealt with Autoregressive Conditional Heteroskedasticity (ARCH) testing in case of subsample changing variances. 33,34 We endeavor in the present context to follow a richer avenue, given the plethora of ad hoc tests and statistics available for deeply scrutinizing the performance of the signal in terms of changing trends, means and/or variances, and cyclical shapes. Among these we test, in a two-sample setting, for major changes that might have occurred within. Specifically, we test from Equation (2) the evolution of the slow-mode term * X t , which includes a linear trend X t  estimable by OLS, and the smoother ˆ* , X t K of Equation (7). We also test for equal variances by means of the F-or Ansari-Bradley test of Equations (3) and (4) respectively. We finally add a test for different subsample shape in terms of cyclical amplitude (e.g. Schwabe cycle).
Let the following: Where the computed TS are penalized whenever one or more of the coefficients j Ψ exceeds or falls short of unity. Penalization consists of the reduction in the absolute value of the computed TS, possibly below the associated critical value or even turning into a positive value. In practice, the farther away are the coefficients j Ψ from their expected values, the greater the risk of nonstationarity in spite of stationarity originally detected by the conventional test statistics. Incidentally, the asymptotic 10%, 5%, and 1% critical values of the computed TS are -2.22, -2.93, and -3.40. [35][36][37] Incidentally, such values are close to some simulated experiments over 10,000 replications of the nonlinear version of Equation (2).

Conclusion
Five statistical practitioner's procedures in support of signal analysis have been introduced and discussed. Such snippets were shown to be useful for better and efficient implementation of the corresponding procedures and of some select hypothesis tests existing in the received literature. Of notable consideration are the optimal smoothing and the corrected ADF test procedures, both characterized by significant novelty in their implementation techniques. The Bayesian stepwise procedures for Monte Carlo simulation and calibration are also of interest because of their ability to perform back/forecasting of signals over long periods of time, if necessary. Of the two, the former may be preferred because quick and easy to implement. The first and the last snippets were found to be particularly useful in a highly erratic environment regarding the signal. The latter, in particular, is worth of notice because of its capability to reverse the results of the traditional stationarity test statistics.


extfieldsofstudy:
Computer Science
Medicine

text:
Multi-Lag Analysis of Symbolic Entropies on EEG Recordings for Distress Recognition

Distress is a critical problem in developed societies given its long-term negative effects on physical and mental health. The interest in studying this emotion has notably increased during last years, being electroencephalography (EEG) signals preferred over other physiological variables in this research field. In addition, the non-stationary nature of brain dynamics has impulsed the use of non-linear metrics, such as symbolic entropies in brain signal analysis. Thus, the influence of time-lag on brain patterns assessment has not been tested. Hence, in the present study two permutation entropies denominated Delayed Permutation Entropy and Permutation Min-Entropy have been computed for the first time at different time-lags to discern between emotional states of calmness and distress from EEG signals. Moreover, a number of curve-related features were also calculated to assess brain dynamics across different temporal intervals. Complementary information among these variables was studied through sequential forward selection and 10-fold cross-validation approaches. According to the results obtained, the multi-lag entropy analysis has been able to reveal new significant insights so far undiscovered, thus notably improving the process of distress recognition from EEG recordings.

INTRODUCTION
Emotions are essential in human communication and interaction, and considerably influence on daily tasks related to cognition, perception and rational decision-making processes (Coan and Allen, 2007). Traditional techniques for emotion recognition are mainly focused on the analysis of physical aspects like facial expressions and speech characteristics (Calvo and D'Mello, 2010). However, given that emotional responses are initiated in the brain and then spread to other biological systems (Gao et al., 2015), interest in electroencephalogram (EEG) signals for emotion recognition has notably increased during the last years Fernández-Sotos et al., 2018;Ramirez et al., 2018).
Existing affect models include from a few basic emotions (Ekman, 1992) to a wide variety of emotional states derived from the combination of basic ones (Schröder and Cowie, 2011). Russell's circumplex affect model is one of the approaches most widely used for emotion classification (Russell, 1980). In this bidimensional approach, emotions are distributed according to their level of valence (ranging from negative to positive) and arousal (from deactivated to activated). A relevant emotion that is receiving growing attention is negative stress, also called distress, because it presents a high prevalence in developed countries (Bong et al., 2013;Alberdi et al., 2016). Although today short-term distress is not considered a risk factor for health, a chronic condition of this emotion often causes or aggravates physical and mental disorders (Bender and Alloy, 2011;Mozos et al., 2017). In this regard, automatic distress identification from EEG signals would help prevent health problems and improve people's quality of life.
Since neural processes are non-linear and non-stationary, both at cellular and global level , nonlinear metrics applied to EEG signal analysis should provide more relevant findings than linear indices traditionally used (Valenza et al., 2012). But, few studies have applied non-linearity to automatic detection of negative stress through EEG recordings (García-Martínez et al., 2019a). This is the case of symbolic entropies, such as Permutation Entropy (PE) (Bandt and Pompe, 2002) and Amplitude-Aware Permutation Entropy (AAPE) (Azami and Escudero, 2016), having demonstrated their efficiency in discriminating between calmness and distress (Hosseini and Naghibi-Sistani, 2011;García-Martínez et al., 2017;Martínez-Rodrigo et al., 2019). Here, the quantification of similar patterns is typically obtained through consecutive samples, or their averaging, within a complete time series.
No lag or time delay between patterns is necessary in those cases where the autocorrelation function of the signal presents a steep decay. However, a time series with long-range linear correlations shows a slow decay in its autocorrelation function. Not applying a lag may hinder entropy metrics from properly quantifying the complexity and non-linear dynamics of the signal. Indeed, it has already been demonstrated that timedelayed entropy tests are helpful to diminish the influence of autocorrelation for better evaluation of the non-linearity of time series (Kaffashi et al., 2008). Hence, a multi-lag approach has been applied to localization of epileptogenic areas through EEG recordings (Zhu et al., 2015).
Let us highlight that an improvement of PE called Permutation Min-Entropy (PME) has been recently introduced (Zunino et al., 2015). PME consists of an improved time-delayed symbolic alternative for identifying the existence of hidden temporal correlations in time series. This allows a better discrimination of time series with similar temporal correlations. Moreover, PME has been very recently applied to emotion recognition by using heart rate variability (Xia et al., 2018). The promising outcomes open a door to the hypothesis that time-delayed analysis may uncover existing information in physiological systems, not revealed before through non-delayed or basic multiscale entropy (MSE) analyses. Furthermore, to the best of our knowledge, no previous research has focused on the study of multi-lag approaches for emotion recognition from EEG signals.
For this reason, in the present manuscript a time-delayed version of AAPE-called Delayed Permutation Entropy (DPE)and PME metrics are applied for the first time with several time delays for the sake of checking the influence of the lag on discrimination between calmness and distress from EEG recordings. The remainder of this paper is structured as follows. Section 2 details the analyzed database, the DPE and PME metrics computed from the EEG recordings and the statistical analysis. Section 3 summarizes the results, which are then discussed in section 4. Finally, section 5 concludes the most remarkable findings related to this study.

Database
EEG signals were extracted from the publicly available Database for Emotion Analysis using Physiological Signals (DEAP) (Koelstra et al., 2012) to guarantee the reproducibility of this study as well as its fair comparison with previous or future works. This dataset contains a total of 1,280 EEG recordings and other peripheral variables from 32 healthy participants with ages ranging 19-37 (mean age of 26.9; 50% male) under different affective conditions. Forty one-minute length video clips with emotional content were used as stimuli in the experiment leading to the dataset. After each visualization, the participants described their emotional state by means of selfassessment manikins (SAM), two graphical scales representing nine intensity levels of valence and arousal (Morris, 1995).
Although the trials contained within the dataset cover the whole valence-arousal space, only two subsets corresponding to distress and calmness emotional states were studied in the present study, as shown in Figure 1. Indeed, calmness and distress groups were selected according to previous works dealing with the same problem (Bastos Filho et al., 2012;Pomer-Escher et al., 2014;García-Martínez et al., 2016;García-Martínez et al., 2017). Hence, distress trials were selected from arousal and valence levels higher than 5 and lower than 3, respectively. On the other hand, the calmness group contained trials with arousal and valence values lower than 4 and between 4 and 6, respectively. Therefore, a total number of 122 and 137 trials of distress and calmness, respectively, were finally analyzed in this work. Moreover, it is important to highlight that only the last 30 s of each trial were selected for further analysis.

EEG Signal Preprocessing
EEG signals were recorded at a sampling rate of 512 Hz with 32 electrodes placed according to the 10-20 standard system of electrode location (Klem et al., 1999). Before starting any kind of analysis, the recordings were preprocessed to eliminate noise and artifacts, thus preserving only the information related to brain activity. To this respect, the signals were initially downsampled to 128 Hz and all EEG channels were re-referenced to the average potential of all electrodes. Next, a forward/backward high-pass filter at 3 Hz and a low-pass filter at 45 Hz cutoff frequency were applied to remove baseline and power line interferences, maintaining the frequency bands of interest in EEG recordings (Koelstra et al., 2012). After that, artifacts derived from physical activity (e.g., facial movements, eye blinks, heart bumping, etc.) and technical sources, such as electrode-pops were eliminated by means of independent component analysis (ICA) (Goh et al., 2017).
A well-known method with ability to automatically identify noisy independent components (ICs) was used (Nolan et al., 2010). Briefly, the algorithm firstly computed correlation between all ICs and electrooculography channels, as well as spatial kurtosis, power spectrum slope, Hurst exponent and median gradient for all ICs. Those components presenting at least an index with a value three times higher than standard deviation for all ICs were then removed. As a final step, the denoised EEG signal was reconstructed from the remaining ICs. It is worth noting that 1.05 ± 0.60 ICs were removed in average for each trial. More precisely, any artifactual IC was identified in 38 trials (14.67%), only one was removed in 168 trials (64.86%) and two were rejected in the remaining 53 trials (20.47%).

Time-Delayed Version of Amplitude-Aware Permutation Entropy
Amplitude-Aware Permutation Entropy (AAPE) is an improvement of Permutation Entropy (PE) to consider amplitude information from analyzed time series (Fadlallah et al., 2013). Although this index has been mainly used in single-lag analyses, it can be adapted to deal with different time scales by changing the embedding delay τ (Azami and Escudero, 2016). Thus, for delayed-time PE computation, a time series Each vector X τ i,m is associated with an ordinal pattern, described as permutation κ i = {r 0 , r 1 , . . . , r m−1 } of {0, 1, . . . , m − 1}, such that its single components fulfill Hence, a total number of m! ordinal sequences π k are obtained from patterns X τ m . Then, the relative frequency of each sequence π k is used to estimate its probability of appearance such that being δ(u, v) the Kronecker delta function modified specifically to work with sequences, i.e., δ(u, v) = 1, if u(i) = v(i), for every i = 1, 2, . . . , m; 0, for otherwise.
Then, delayed-time PE is finally obtained by computing the Shannon entropy from the probability distribution of all symbols, such that The index is normalized by term ln(m!) to obtain values ranging between 0 and 1. In the case of a completely predictable signal, only a pattern π k is found and PE reports a 0 value. On the contrary, symbols π k in unpredictable time series present the same probability of occurrence. Thus, PE provides the highest value 1. Hence, predictability information reported by PE is easily interpretable (Zanin et al., 2012). Nevertheless, only the ordinal structure of patterns is considered by this index, thus discarding the information related to the amplitude of each sample. As amplitude differences could play a key role to determine the predictability of a time series, AAPE was introduced to overcome this limitation (Azami and Escudero, 2016). AAPE computation is based on calculating the probability of repetition of each pattern π k by considering its relative frequency, and also the average absolute (AA) and relative amplitudes (RA) of vectors X τ i,m . Amplitudes AA and RA are obtained, respectively, for a specific vector X τ i,m as and Then, the relative frequency of π k is computed as being K an adjusting coefficient of terms AA and RA, ranging from 0 to 1. As recommended by the authors, a value K = 0.5 was considered here. Finally, delayed-time AAPE, referred to as Delayed Permutation Entropy (DPE), is computed by means of Shannon entropy, such that

Permutation Min-Entropy
Recently, PE has also been generalized by replacing Shannon entropy with Rényi one, reaching a better characterization of some rare and frequent ordinal patterns . More precisely, Rényi Permutation Entropy (RPE) is defined as where order q (q ≥ 0 and q = 1) is a bias parameter. Indeed, q < 1 benefits rare events and, contrarily, q > 1 privileges salient ones. It is mandatory to note that Shannon entropy is an instance of Rényi entropy for q = 1 and, hence, RPE is a more flexible tool than PE. In this respect, RPE has reported a more complete characterization of a variety of complex dynamics, including physiological processes (Mammone et al., 2015). In addition, RPE is featured to converge to a minimum entropy in the limit q → ∞, thus providing Permutation Min-Entropy (PME) (Zunino et al., 2015). This new entropy-based metric is quickly and simply computed as still retaining the main advantages of PE, such as its simplicity, low computational cost, noise robustness, and invariance with respect to non-linear monotonous transformations. This index has also proven a greater ability than PE to detect the existence of subtle temporal structures in EEG channels (Zunino et al., 2015).

Feature Selection
Firstly, single DPE and PME values from lag τ = 1 to lag τ = 10 were computed for each subject by using a pattern length m = 6. Entropies computed for each time-lag are a measure of predictability of the time series and assess the effect of time dynamics from an inter-lag perspective. Indeed, larger entropy values represent more unpredictable dynamics of the EEG signals, showing an increase of autocorrelated patterns in a long-term fashion.
On the other hand, multi-lag entropy curves were parameterized by means of slopes, areas under curves and arc lengths. Indeed, some studies have previously reported that features extracted from parameterized curves may reveal important information related to the dynamics of the signals across different temporal intervals (Escudero et al., 2006). In this regard, to estimate the trend evolution of each time-lag curve, slopes between delay τ = 1 and τ = 2, 4, 6, 8, and 10 were calculated from all EEG channels of each trial and denoted as Slp1 − τ . The slope is estimated as the straight line connecting the multi-lag entropy values under study. Higher slope values suggest larger entropy increases between the original signal (τ = 1) and higher versions in consecutive multi-lag time delays (τ = 2-10).
Furthermore, areas enclosed under the multi-lag curve between lag τ = 1 and lags τ = 2, τ = 4, τ = 6, τ = 8, and τ = 10 were computed and denoted as Ar1 − τ . In this sense, a higher area is achieved when DPE and PME entropy values are higher for the majority of time delays, suggesting that time series are less predictable. Finally, the arc length (AL) for each time-delayed curve was computed between lags τ = 2 and τ = 10. An arc length value shows the morphological alterations of the curve across different lags, and may show significant differences among lags from different groups of study. The arc length of each multi-lag curve was computed as referring E to the values of DPE and PME for the corresponding time-lag τ in each case. Hence, a total of 21 features were computed for symbolicbased entropies DPE and PME on each EEG channel. More precisely, 10 single entropy values (one for each of the 10 timelags computed), 5 tendency parameters related to time-lag curves (slopes Slp1 − 2, Slp1 − 4, Slp1 − 6, Slp1 − 8, Slp1 − 10), and 6 shape-related features (areas under curves Ar1 − 2, Ar1 − 4, Ar1 − 6, Ar1 − 8, Ar1 − 10, and arc length AL) were obtained for each EEG channel.

Statistical Analysis
Once the features were computed for each metric under study, Shapiro-Wilks and Levene tests corroborated the normality and Frontiers in Neuroinformatics | www.frontiersin.org homoscedasticity of the data, such that the results are expressed as mean and standard deviation. Then, statistical differences between features obtained for emotional states of calmness and distress were assessed for each time lag τ using a one-way analysis of variance (ANOVA). A value of statistical significance ρ < 0.05 was considered as significant.
Furthermore, the discriminatory power of each feature to distinguish between both groups of emotions was tested by using a stratified 10-fold cross-validation scheme. This methodology prevents over-fitting as well as other biases when performing the training/test operation on classifiers (Jung and Hu, 2015). Thus, the database selection containing 259 recordings was sliced into ten equally-sized folds with a balanced number of trials from both groups. Next, ten iterations were performed, such that in each one 9 out of 10-folds were used as a training subset, and 1 out of 10-folds was used as the test subset. To perform the classification, a receiver operating characteristic (ROC) approach was computed using the training trials to obtain an optimal threshold, which was then used to classify the trials in the test subset. It is worth noting that the threshold was selected as the cut-off point that maximizes accuracy (Acc). Values of sensitivity (Se), specificity (Sp), and Acc, obtained from the 10 iterations, were finally averaged to provide global and robust estimates.  Keeping in mind the objective of assessing possible relationships and complementary information among features, several advanced classifiers were used. Thus, a decision tree classifier (DTC), a support vector machine (SVM), a quadratic discriminant analysis (QDA) and a k-nearest neighbor (KNN) classifier were used. Regarding DTC, the nodes' growth was stopped when each node solely contained either fragments from only one group or a number of trials <20% of the entire dataset. Moreover, every node was split by using an impurity-based Gini index. Furthermore, the SVM classifier was run with a cubic kernel function and kernel scale of 1. Finally, the KNN classifier used an Euclidean distance metric with 10 neighbors, where the weight of the distance was computed to perform the classification by means of squared inverse. Nonetheless, given the high amount of analyzed features (21 features × 2 metrics × 32 channels), the subset providing most information was selected in first place for each classifier. Thus, a sequential forward selection (SFS) approach was used to select the subset of features minimizing misclassification rate for each classifier. A stratified 10-fold cross-validation scheme was also used to reduce overfitting in this analysis.

Results for Delayed Permutation Entropy and Permutation Min-Entropy
Mean and standard deviation of DPE and PME values for the most relevant EEG channels at different time-lags (1 ≤ τ ≤ 10) are shown in Figures 2, 3, respectively. As can be observed, both metrics obtained a similar trend throughout the increasing time-lags. DPE and PME values for calmness are higher than for distress trials, especially at lower lags. However, as time-lag increases the average differences between groups become smaller, such that at higher time-lags the mean entropy differences between groups become imperceptible. Furthermore, a certain degree of stabilization at time-lags > 3 for both metrics can also be noticed, where the standard deviation decreases as the analyzed time-lag increases.  Table 1 shows the statistical significance and global classification performance, that is accuracy (Acc), for timelags (1 ≤ τ ≤ 10) in DPE and PME metrics, respectively. Although there are similarities among mean entropy curves, the performance achieved for each metric differs considerably throughout the time-lags. In general terms, DPE shows a poorer performance discriminating between emotional states of calmness and distress than PME. As can be observed in Table 1, only lower time-lag entropies show a relevant statistical significance. Moreover, the overall discriminatory power for all cases is around 60%, decreasing even more when higher time-lags are analyzed. This effect is clearly seen at parietooccipital and occipital channels PO3, O1, O2, and Oz. For instance, PO3 channel achieved a global performance of 62.9% at time-lag τ = 1, while accuracy decreased down to 57.10% at time-lag τ = 9. Only parietal channel P3 showed a regular statistical significance throughout every single time-lag, reaching a maximum global classification performance at time-lag τ = 2 with 63.7% of subjects classified correctly.
On the contrary, PME showed better global classification performance for certain channels, especially at higher time-lags. For instance, parietal channel P3, classified correctly 68.3% of trials between calmness and distress at time-lag τ = 9. Similarly, centro-parietal channel CP2 showed a poor performance when no lag was applied (τ = 1), but it raised at higher time-lags, achieving the best single global performance at time-lag τ = 9, classifying correctly the 68.7% of the subjects. This improvement supposes an increase of more than 13% regarding PME at no lag τ = 1 and more than 8% compared with DPE metric at time-lag τ = 9. Moreover, there seems to be a certain degree of complementarity between DPE and PME at different time-lags, because the same EEG channels measured with each metric show relevant information at different time-lags. This contrast is wellnoticed at parietal channel P3. No relevant differences between several time-lags were found when DPE was computed for P3, i.e., all time-lags presented a similar discriminatory power. On the contrary, the same channel showed an important increasing performance when it was analyzed by means of PME metric at higher time-lags; hence relevant information was noticed when time-delay was performed. Table 2 summarizes diagnostic accuracy of every curve-based parameter derived from DPE analysis for the most relevant channels. As can be seen, almost all features achieved statistically significant differences between groups (ρ < 0.05). Features obtained from parietal channel P3 achieved a notable statistical significance, especially in the area under the time-lag curve (Ar1 − 2 to Ar1 − 6), thus stating the differences between curves at lower and their convergence at higher time-lags. It is also remarkable that parieto-occipital channel PO3 achieved a good performance for all features. With respect to slopebased parameters, global accuracy ranged from 60.62 to 62.23%, Slp1 − 10 being the feature with maximum performance for this channel. Similarly, accuracy of area-based lag parameters ranged from 59.07 to 62.16%, where Ar1 − 2 reported the maximal performance. It is also worth noting that arc-length reached the maximum global accuracy, classifying correctly 63.47% of trials, and thus overcoming the best performance obtained by single DPE entropy at the same channel at lag τ = 1. Finally, the rest of the parameters computed from occipital channels O1, O2, and Oz and centro-parietal channel CP2 obtained a more limited performance, their global accuracy ranging from 53.28 to 61.12%. Similarly, Table 3 summarizes discriminant ability of all curve-based parameters derived from PME analysis for the most relevant channels. In this case, both statistical significance and global accuracy are more limited than for DPE curve-related parameters. All features computed on parietal channel P3 and parieto-occipital channel PO3 resulted to be statistical significant. The global accuracy obtained for these channels ranged from 57 to 61%, thus achieving a worse performance than before. Moreover, only a few curve-related parameters from occipital channels O1, O2, and Oz showed statistical significance, and the global accuracy was below 60% for all the parameters. Finally, CP2 achieved the worst performance, where global accuracy was around 55%.

Multi-Parametric Analysis and Advanced Classification
For each classifier, the optimal number of features minimizing its misclassification rate through an SFS scheme ranged from 5 to 8 in each iteration of a 10-fold cross-validation approach. The occurrence of the most relevant variables are displayed in Figure 4. As can be seen, entropy-based metrics were mainly  chosen for time-lags longer than 1 and curve-related variables both from DPE and PME. Moreover, it should be noted that most of these features were selected in nearly all folds, thus only changing the less relevant ones for the resulting classification models. More precisely, for all classifiers, most DPE-based parameters were chosen from EEG channels P3 and PO3 and most PME-based features from channels P3, Pz, FC5, C4, and CP5. Once the feature selection process finished for each classifier, the obtained classification results are displayed in Table 4 in terms of sensitivity, specificity and accuracy. Note that global discriminant ability for all approaches ranged from 75.66% (for QDA) to 92.32% (for KNN). Furthermore, the SVM classifier achieved a comparable performance to KNN, classifying correctly 91.12% trials. Finally, it should also be highlighted that both SVM and KNN classifiers reported the largest diagnostic accuracies with well-balanced values of sensitivity and specificity.

DISCUSSION
During the last years automatic emotion recognition has received special attention due to its importance in areas, such as medicine and education. Among the different types of emotions, continued distress is one of the most studied because it is often harmful for health. Considering its relevance, distress has been assessed in a wide variety of scenarios, including driving tasks (Healey and Picard, 2005), military exercises (Skinner and Simpson, 2002), surgical procedures (Marrelli et al., 2014), and on-line exams (Gomes et al., 2014), among others. An interesting study recently published shows a methodology to redirect stress episodes toward positive moods (Fernández-Caballero et al., 2016).
Taking into consideration this preamble, several works have been published in the literature. Their research is focused on FIGURE 5 | Topological EEG representation of average calmness and distress values for time-lags from τ = 1 to τ = 9 using PME metric.
automatic distress recognition using EEG recordings (Hosseini et al., 2010;Khosrowabadi et al., 2011;Peng et al., 2013;Minguillon et al., 2016;Al-Nafjan et al., 2017;Al-Shargie et al., 2018;Jebelli et al., 2018;García-Martínez et al., 2019b). However, only a few of them have analyzed this phenomena from a non-linear point of view (García-Martínez et al., 2016). Recently, another approach reported that symbolic analysis of brain dynamics was able to detect distress (García-Martínez et al., 2017). In that study, PE and its extension called AAPE were used to assess brain dynamics for each EEG channel. However, the analysis was carried out without considering the possibility of exploring hidden non-linear information at timelags higher than one. This was the starting point that motivated the present study.
To the best of our knowledge, this is the first work addressing the effects of multi-lag for distress recognition from EEG recordings. For this purpose, a modified version of AAPE was used to analyze EEG signals with distinct timelags. Additionally, PME was also considered in this study, since it is an improved symbolic alternative for identifying the existence of hidden temporal dynamics in time series and it allows a better discrimination of signals with similar temporal correlations (Zunino et al., 2015). Indeed, it has been recently applied in the study of emotion recognition using heart rate variability with promising results (Xia et al., 2018).
As expected, both DPE and PME metrics reported the same trends when calculating the mean entropy values across the different time-lags, as was observed in Figures 2, 3. Calmness emotional state reported higher entropy values than distress for all EEG channels, especially at lower time-lags. However, this difference became smaller as the time-lag increased. Furthermore, the time-lag analysis revealed additional entropy information not observed at time-lag τ = 1, especially at centroparietal and occipital channels. This effect can also be wellnoticed in Figure 5, which shows a topological representation of mean entropy values for each channel at the first nine timelags using PME. Although the general trend is maintained across time-lags, the imprint patterns change, thus revealing information at certain channels not seen before.
These results enhance the presence of a larger diversity of ordinal patterns in some local time series in calmness trials in comparison to distressed ones, thus suggesting the existence of more complex brain dynamics in calmness state. Such loss of complexity under distress might be associated to a lower brain's ability of adaptation to external stimuli and environmental changes. Indeed, decomplexification of physiological systems has been traditionally identified with a lower ability to manage information, and therefore with a higher probability of suffering a pathological condition (Goldberger et al., 2002;Lipsitz, 2004). Interestingly, these findings are in agreement with other studies FIGURE 6 | Topological EEG representation of average calm and distressed patients for time-lag τ = 1 vs. τ = 9 using PME metric.
published during the past years. In this sense, increased values of correlation dimension in calm participants with respect to distressed subjects have been reported so far (Hosseini et al., 2015). In another work, a decrease of relative power in subjects facing a series of distressful stimuli was also described (Bastos Filho et al., 2012). Finally, the fact that this trend is maintained throughout the different time-lags reinforces a previous study where no lag was applied (García-Martínez et al., 2017).
Another relevant finding is that single discriminatory power in multi-lag analysis has notably improved with respect to other previous studies dealing with singe-lag (García-Martínez et al., 2017) and MSE (Martínez-Rodrigo et al., 2019) analysis, especially in some specific channels. Thus, left parietal channel P3 is still a very relevant channel for distress detection using symbolic analysis. This finding was already reported in our recent previous work where AAPE was applied to the data (García-Martínez et al., 2017). Nevertheless, other studies have already corroborated this association with the left parietal area. Thus, a higher activation has been observed during normal nondepressed and reasonably positive moods in the left parietal area than in the right one (Davidson, 1988). In the same line, meditation has also been characterized by an increasing activity of the left parietal region (Manna et al., 2010). In the present study, P3 showed robustness and consistency across different time-lags when discriminating between emotional states of calmness and distress. Nevertheless, the global classification performance was improved notably for time-lags higher than one, especially when data was analyzed by means of PME metric, increasing from 61% when no lag was applied up to 68.30% when time-lag τ = 9 was computed.
These findings may indicate the existence of long-range correlations in the data, which have only been sufficiently highlighted by considering a multi-lag entropy-based analysis. Indeed, these observations can be visually corroborated in the topological representation of brain areas depicted in Figure 6. It represents average PME values computed from emotional states of calmness and distress for time-lags τ = 1 (a) and τ = 9 (b). As can be seen, entropy values obtained at time-lag τ = 9 are quite different compared to the analysis with no lag, showing a more balanced pattern between left and right hemispheres throughout frontal, parietal and occipital areas. In addition, the entropy differences between calmness and distress are also shown in this figure (on the right column). Thus, the higher differences are found in left central region for τ = 1, whereas a higher activation of left parietal region is obtained for τ = 9.
The right frontal channel F4 also presents a considerable difference of activation between calmness and distress states both in τ = 1 and τ = 9 cases. Interestingly, the relevance of the mentioned areas and the possible relation between frontal and parietal areas of opposite hemispheres has already been depicted in our previous studies. In fact, 30 years ago it was verified that a relative left parietal brain activation is balanced by a relative right parietal brain activation and vice-versa (Davidson, 1988). A similar outcome has also been observed in another study where patients with different mental disorders were conducted to practice meditation (Rubia, 2009).
Moreover, the right brain hemisphere deserves especial attention in this work. Considering our previous findings, the right centro-parietal channel CP2 showed no relevance when analyzing emotional states with symbolic entropies (García-Martínez et al., 2017). This outcome has been corroborated again in this work, where neither DPE nor PME showed statistical significance at τ = 1 (no lag), and the global classification was below 55% in both cases (see Table 1). However, when analyzing the same brain area at higher time-lags, a notable increase of statistical significance and discriminatory power was observed. It is especially the case for CP2 with PME at time-lag τ = 9 (ρ = 0.0026 and Acc = 68.70%), thus achieving the highest global classification in this study. The relevance of the channel CP2 can also be observed in Figure 6, where the difference is not notable for τ = 1, but it is for τ = 9.
These findings, together with the relevance of the results reported by the left parietal channel P3, reveal the possible existence of complementary information among the parietal lobes of both brain hemispheres. Indeed, a number of previous works reported a certain degree of complementarity between right and left posterior areas under stimulation of distress and calmness. For instance, interesting information about parietal and occipital asymmetry at different frequency bands during distressful tasks has been described (Park et al., 2011). Furthermore, an intensive parietal lobe activation under anxiety and distress conditions has also been reported (Nitschke, 1998).
Recently, occipital electrodes O1 and O2 have also been explored to evaluate variations in complexity provoked by visual elicitation (Tonoyan et al., 2016). In a similar line, the posterior brain area has been related to the arousal component of emotions, thus being their processing essential for the recognition of emotions (Dolcos and Cabeza, 2002). Interestingly, in our previous study the combination of the left parietal channel P3 and the right parietal channel P4 achieved a notable performance discriminating between emotional states of calmness and distress, thus demonstrating the inter-correlation of these brain regions (García-Martínez et al., 2017). However, in that study, brain dynamics were assessed by means of different computation approaches, where each methodology highlighted one of the hemispheres in isolation. In the present work, both areas have resulted to be significant when analyzed with PME metric at higher time-lags, obtaining results comparable to those reported individually by other metrics used in our previous works.
The obtained multi-lag curves were also parameterized and studied to compare the relative complexity of normalized time series. The use of curve profiles for characterization of biological signals has been already proposed by other authors (Costa et al.,  2005). Accordingly, Slp1−4 and Slp1−6 reflect that the degree of change in the complexity of some EEG channels is more relevant in smaller time-lags. These outcomes were observed in Figures 2,  3, where changes in the slopes could be seen until the curves stabilized around time-lag τ = 5. In this regard, parietal and occipital channels showed statistical relevance at these slopes for DPE metric, but only parietal channels resulted relevant when curved-related parameters were calculated for PME metric. It is worth noting that the same outcomes were also obtained for the area under curve parameters.
On the other hand, the developed multivariate analysis has shown that putting all the data together led to a notably overall performance increase, which demonstrated that multi-lag analysis is able to provide additional, as well as complementary information, to single-lag one. To this respect, an SFS scheme was applied under a 10-fold cross-validation approach to choose the optimal subset of features maximizing the classification rate for each classifier. Interestingly, the most relevant features selected for each classifier were mainly computed from EEG channels showing the largest differences between emotional states of calmness and distress. Indeed, DPE computed from channels PO3 and P3 reported high statistically significant differences between the two groups of trials, as shown in Tables 1, 2. Similarly, PME obtained from channels P3, Pz, FC5, C4, and CP5 provided high visual differences, as shown in Figures 5,  6. Moreover, let us highlight that all selected features were computed from time-lag longer than 1 or from τ -based curves, as shown in Figure 4.
Another relevant aspect is that the results obtained in the present study outperformed notably other similar works that have analyzed non-linear metrics from the same database with no lag, such as summarized in Table 5. Indeed, a global accuracy of 69.6% has only been reported by applying a high-order crossing approach to four EEG channels (Bastos Filho et al., 2012). In addition, different non-linear metrics have reported a higher level of complexity in stressed subjects (Peng et al., 2013). In other work, combining quadratic sample entropy values from several EEG channels through a DTC classifier, a discriminant ability around 75% has been provided (García-Martínez et al., 2016). On the other hand, a discriminant model based on SVM and using irregularity and symbolic metrics reached a diagnostic accuracy >80% (García-Martínez et al., 2019b). Similarly, variants of PE have already been applied to distress recognition, with a classification performance of 81.31% (García-Martínez et al., 2017). However, in the present study an improvement of about 10% has been reported by making use of the same kind of SVM classifier.
Likewise, the classification results obtained in the present study also significantly improved the sole recent work conducting a MSE analysis on EEG signals for distress identification (Martínez-Rodrigo et al., 2019). In fact, making use of the same SVM-based approach, here a classification rate between distress and calmness emotional states has been obtained about 6% higher than for MSE. It should be noted that both MSE and multi-lag entropy analyses pursue the same goal of quantifying complexity at different time scales. For this purpose, MSE uses a rescaling procedure based on filtering out the shorter oscillations and keeping the longer ones (Humeau-Heurtier, 2015). This approach unavoidably removes some frequency content, specially from rescaled time series presenting very fast oscillations (Humeau-Heurtier, 2015). Such loss of frequency information could explain the aforementioned poorer outcome reached by MSE, because entropy computation from time-lagged samples does not alter time and frequency information from original data (Govindan et al., 2007;Kaffashi et al., 2008). Moreover, this finding could also justify the fact that, whereas no changes were noticed across all time scales in MSE analysis in brain areas activating and supporting distress (Martínez-Rodrigo et al., 2019), large differences have been observed for different time-lags, as extensively described in previous paragraphs.
Finally, there are some limitations in this study that deserve our attention. Firstly, the studied DEAP database is not exclusively designed for recognition of calmness and distress emotions. In fact, many other emotional states were also recorded during the experiment (Koelstra et al., 2012). Moreover, the number of trials eliciting calmness and distress is notably unbalanced for each healthy volunteer, thus making the use of a subject-based classification impossible. Secondly, further analyses on other similar databases like ASCERTAIN (Subramanian et al., 2018), AMIGOS (Miranda Correa et al., 2018), and DREAMER (Katsigiannis and Ramzan, 2018) are required to corroborate and generalize the obtained results. However, the impact of some potential confounding aspects on the results provided by several databases will have to be carefully analyzed for this purpose. Thus, it should be thoughtfully scrutinized how different experimental setups, population distributions in terms of age and gender, and technical aspects related to the acquisition of EEG signals mask changes in brain dynamics under distress. Thirdly, the video clips used as stimuli have a duration of 1 min, which may be too much time to just elicit a single emotional state. Thus, participants may present different emotions for the same stimulus, making it difficult to properly rate their level of valence and arousal. Finally, only EEG signals were assessed in this work, thus discarding the information reported by other physiological variables. However, peripheral recordings also contained in DEAP and other databases, in combination with brain dynamics, will be explored in further studies for the sake of detecting distress episodes.

CONCLUSIONS
In this study, two permutation entropies, adapted to work in a multi-lag context, have been analyzed for the first time to automatically identify negative stress. This multi-lag analysis has revealed new insights never seen before, thus notably improving the performance of distress identification. Considering the relevant results that permutation entropy has previously reported in non-lag and multiscale contexts for human emotion detection, it becomes highly interesting to analyze brain dynamics from a time delay viewpoint. For this reason, we hypothesized that there might exist relevant and complementary information at higher time-lags among different brain areas. The results obtained after performing the analyses have confirmed our initial ideas, reporting an improved classification between emotional states of calmness and distress.
Moreover, left parietal and right centro-parietal channels showed remarkable activation at higher time lags, suggesting that removing long-range linear correlations may help to better evaluate the non-linear information of the data. Finally, several discriminant models obtained from advanced classifiers were used to study the complementarity of the features computed at different time-lags for each EEG channel. The resulting functions have combined single entropy values from different channels calculated at lags higher than one with curve-related parameters, thus corroborating that there are more relevant information when time-lags are applied to the time series than when data are analyzed without any time delay or averaging consecutive samples.

AUTHOR CONTRIBUTIONS
AM-R and BG-M conceived and designed the study, programmed the experiments and drafted the manuscript. LZ and RA helped to interpret the results and reviewed the manuscript. Finally, AF-C supervised the experiments, reviewed the manuscript and contributed to the final version. All authors have read and approved the final version of the manuscript.


A Novel Online Incremental Learning Intrusion Prevention System

Attack vectors are continuously evolving in order to evade Intrusion Detection systems. Internet of Things (IoT) environments, while beneficial for the IT ecosystem, suffer from inherent hardware limitations, which restrict their ability to implement comprehensive security measures and increase their exposure to vulnerability attacks. This paper proposes a novel Network Intrusion Prevention System that utilises a SelfOrganizing Incremental Neural Network along with a Support Vector Machine. Due to its structure, the proposed system provides a security solution that does not rely on signatures or rules and is capable to mitigate known and unknown attacks in real-time with high accuracy. Based on our experimental results with the NSL KDD dataset, the proposed framework can achieve on-line updated incremental learning, making it suitable for efficient and scalable industrial applications.

I. INTRODUCTION
Commercial intrusion detection/prevention systems typically suffer from low detection rates and high false positives which require substantial optimization and network specific fine tuning. The majority of these systems rely on signatures to detect potential attacks and, therefore, cannot detect "zero day", unknown attacks. Despite the significant research effort to introduce intelligence to NIDS/NIPS by implementing anomaly based/machine learning detection methods, there has been little real-world adoption of such systems. A machine learning-based NIDS/NIPS partially addresses the above inherent limitation that signature-and rule-based industrial NIDS/NIPS suffer from. Beyond the lack of real-world machine learning-based applications, there has been very limited academic work focusing on incremental learning algorithms applied to NIDS/NIPS. Incremental learning algorithms allow the classifier to refine and improve its capabilities over time (as it processes increasing amounts of input data) in contrast to an offline or batch learning algorithm, where the classifier is assumed to be exposed to the input data in a batch. Network data dynamics change significantly over time and applying static learned models progressively degrades the detection performance, making offline algorithms unsuitable for a NIDS/NIPS. We propose a novel Network Intrusion Prevention System that exploits the benefits of incremental machine learning frameworks and achieves accurate results, comparable to most offline neural network-based NIDS. The framework is based on a modified version of Self Organizing Incremental Neural Network (SOINN) [1] to achieve on-line clustering, coupled with multiple SVMs to perform classification. The evaluation of the resulting framework was performed on the NSL-KDD dataset [2], which is an improved version of well-known KDD'99 dataset. In spite of its age, the dataset is still the defacto alternative for benchmarking techniques and tools that aim to provide effective intrusion detection systems [2]. Given its widespread use, which makes it easier to provide a reference analysis, we have also adopted this dataset to initially test our method. The results show that the proposed framework can achieve on-line updated incremental learning in a fast and efficient manner. The rest of the paper is organized as follows: the related work is discussed in Section 2, then Section 3 presents the proposed method. Section 4 gives an overview of the experimental results. Section 5 discusses the experimental results and Section 6 the provides a set of conclusions and future work to extend the proposed framework.

II. RELATED WORKS
Intrusion Detection Systems can be categorized based on the detection method used. According to Liao et al [3], detection methodologies are classified in two categories: signaturebased and anomaly-based. The signature-based detection methodology is using patterns previously defined by known attacks and is distributed as a set of signatures. The signatures are then compared against patterns found in the network traffic to discover possible attacks. While effective for known threats, such a method cannot discover or prevent unknown attacks and is unable to maintain and update signatures for known or newly discovered attacks. In contrast, anomalybased detection typically establishes a baseline/normal level using statistically significant traffic. Depending on the data analysis technique used, the training/testing process can use either classification or clustering. There has been substantial work and research towards improving the classification techniques for Intrusion Detection Systems, mostly focusing on evaluating alternative solutions for the core analysis, including neural networks [4] [5], fuzzy logic [5] [6], genetic algorithms [8], and support vector machines [9]. At the same time, clustering techniques mainly exploited the k-means and outlier detection algorithms [10], [11]. As explored in our proposed approach, prior research demonstrated the benefits of combining both stages in a hybrid approach. In their papers [2] [12] combined anomaly and misuse detection to propose a novel solution for an IDS. Following a different combination, the authors in [13] extended a neural based technique with SVMs to propose a solution for an IDS that could recognize both anomalies and known intrusions. Artificial Neural Networks were also evaluated as a possible analysis tool for Intrusion Detection Systems; in [13], the authors proposed a user behaviour model where each user is represented by a neural network, essentially a complement of a statistical model. Later, in [14], the authors used the back propagation algorithm to train the neural network for an IDS proposal with a 96% detection accuracy and a 7% false alarm;  in [15] built a Hierarchical Network Intrusion Detection system using statistical pre-processing and neural network classification. Following from prior studies, [16] showed that Multilayer Perceptron neural networks used for offline analysis could be applied in intrusion detection not only to recognise attacks but also to classify types of attacks. This approach, using two hidden layers, led to a 91% accuracy. Authors in [17] used a distributed time delay neural network to solve a similar multi class problem but with improved accuracy results of 97.24%. Xiang et al used a modification of SOINN [2] and applied it to intrusion detection [18]. With a semi-supervised learningbased IDS proposal, they demonstrated that the user input could be minimized by combining a modified SOINN and SVM to achieve semi-supervised learning with the same space efficiency as a supervised SVM.
A number of more recent studies extended the SOINN approach for a range of IDS-related tasks. To begin with, in [19], the authors proposed an Improved Incremental Learning SOINN for cloud-based environments. While the work aimed to detect functionality anomalies, it also provided a starting point for this work. SOINN can also be successfully applied to visualisation, as demonstrated in [20]. SOINN was also tested on an older dataset (KDD 99) in [21] and led to relatively good accuracy with only minimal pre-processing and feature extraction. On the implementation side, [22] demonstrated the efficiency of using TensorFlow towards analysing the NSL-KDD dataset and detecting malicious traffic with 96% accuracy.
As highlighted by prior research, SOINN and Incremental Learning are indeed very effective methods to tackle the challenges of intrusion detection. We aim in this paper to extend prior work by using incremental machine learning to overcome the issues from using static (offline) machine learning models applied to Network Intrusion Detection/Prevention Systems designs.

III. PROPOSED METHOD
Throughout the evolution of IDS, it became apparent that a hybrid approach leads to better results due to its ability to discriminate using initially statistical-based approaches to improve the formulation and formatting of inputs, then apply them onto an artificial intelligence-based engine. Based on this approach, this paper proposes a conceptual framework designed to monitor inline the observed traffic and process the raw traffic to extract a set of statistical parameters that can then be passed through a detection engine. Once the detection engine determines whether the traffic is part of an attack, the packets can be logged/dropped or passed to the recipient devices. The following subsections present the overall architecture and encompassing blocks.

A. n-SOINN-WTA-SVM
The framework, as seen in Fig. 1, consists of a detection engine that represents the core of the framework, a preprocessing module for the incoming traffic that feeds the detection engine, a validation module to evaluate the results of the detection engine and an update module that feeds back the failed results to the detection engine At the core of the architecture is the detection engine, named "n-SOINN-WTA-SVM". The pre-processing module captures, extracts, and calculates the TCP-related parameters that are used as inputs for the detection engine. The validation module evaluates (either manually or automatically) the results of the detection engine and saves them in order to feed the detection engine in the update phase, whether correct or incorrect. After the initial training, the system switches between two different phases: the live phase, where the input traffic is categorized as normal or attack, allowing the system to take appropriate action, and the updating phase, where the system is incrementally updating its machine learning capabilities by learning from new input data, more specifically failed classifications. Our dataset included four different categories of attacks fed to the detection engine during experimentation and evaluation.

1) Concept.
The core concept of the proposed framework is to incrementally build a network protection mechanism. In the initial phase, the detection mechanism learns using a relatively small sample of network data, sufficient for basic network protection. Subsequently, as more network data becomes available, the mechanism is incrementally updated with input data of classes that it failed to detect, in order to refine and advance its protective capabilities. To support the learning process, the validation mechanism decides whether a decision failed or not. In order to advance its capabilities, the core detection engine must be able to categorize the network data in a multi-class manner, not just whether a connection is an attack or not, but what kind of attack that is and therefore provide a solution to an incremental learning multi-class problem.

2) SOINN.
The self-organizing and incremental neural network (SOINN) [1] is an on-line unsupervised learning mechanism for unlabelled data. SOINN has already been used in other studies as a clustering method that handles supervised data. As an online incremental clustering method, SOINN offers relatively high computational speed with low computational cost. Furthermore, the network complexity and size of SOINN are controlled and stabilised through a "garbage collector" technique. The technique defines an attrition parameter called age which represents the time after which period nodes will be removed if they are not updated for a specified time, and thus dynamically eliminate noise in the data. This property makes it attractive for dynamic environments where long-term learning is required. To ensure its scalability when expanding, the size and growth of the network are controlled by a parameter n where multiple pairs of SOINNs will be used as a supervised clustering method. SOINN, as summarised in Fig. 2, initializes the network with an empty set of nodes, then adds the first two nodes to the list, with weight vectors set as the two input vectors. After the initialization, the neural network finds for every input vector the nearest node (winner) and the second-nearest node (second winner) of the input vector by measuring the distance S1 and S2 from every input to every node with the equations (1) and (2): If the input vector belongs to the same cluster as the winner or second winner based on the distances calculated against a similarity threshold criterion, SOINN updates the weight vector of the node and its neighbours with the weight vector of the input vector and connects it to the node by an edge. If the input vector does not belong to the same cluster of the winner or second winner, the mechanism adds a new node to the network.
3) n-SOINN. The framework proposed by this paper exploits the concept behind n-SOINN [23], which modified the original SOINN to utilize multiple pairs of SOINNs using a supervised clustering approach. n-SOINN uses two significant modifications: a global parameter that controls the topology of the network and uses the squared (rather than the basic) Euclidean distance to calculate the distance between the input and the nodes. In order to control the number of output nodes of the network, where a difference of how accurate the compressed information will be created, a parameter named n is introduced. This parameter dictates that any first winner node that wins more than n times assign a win to the second winner node. If the second winner node also has more than n winning times, a new node is generated. For n=0, the network behaves exactly like the original SOINN. Setting n to a very high values reduces the number of created, stable nodes and favours only the popular ones, which is likely to lead to poorer accuracy. The Euclidean distance used in the original SOINN was intended for the purpose of a single SOINN to realize the unsupervised incremental learning task. The use of the squared Euclidian distance in n-SOINN allows measurement of the distance between nodes for different SOINNs; the aging built in the SOINNs garbage collector allows eliminating nodes of the input data whenever they become unpopular.

4) WTA-SVM.
A Support Vector Machine (SVM) is a supervised machine learning model that operates as a group classifier by constructing a hyperplane or set of hyperplanes in a high dimensional space using training data to separate data into two groups. One of the main advantages of SVM is its effectiveness in high dimensional spaces and its ability to discriminate data which are not readily separable by simpler methods. If the 'normal' and 'abnormal' data sets to be discriminated are not linearly separable or the variations of features for two classes overlap, the SVM maps the original data space into much higher space to make the separation easier. The function for mapping of features and the other parameters should be optimized with different optimization methods. An SVM acts as a binary classification algorithm that solves binary problems. In order to solve a multi-class problem with SVMs, the common approach is to reduce the single multiclass problem to a set of multiple binary classification problems. Two common methods are employed to solve the problem: the one-versus-all method, using the winner-takesall strategy, in which the classifier with the highest output function gets the class assignment, and the one-versus-one method using a max-wins voting strategy, in which the class with the highest number of votes in every win of the two class problem gets the class assignment. This paper uses the former approach (one-versus-all) method and the winner-takes-all (WTA) strategy is employed to solve the multi-class problem.

5) n-SOINN-WTA-SVM.
The detection engine and the core of the framework, as seen in Fig. 3, is initialized with a dataset of k attack classes where x is a d-dimensional TCP feature vector of the connection and y is an attack class category label. A network connection attack category is modelled by two n-SOINNs, one with a high n value and another with a low n value. For each class k, a pair of two n-SOINNs and one SVM binary classifier are created. The n-SOINN with the low n value is supposed to hold more accurate compressed information than the n-SOINN with the high n value, making it a binary problem.
For every pair of n-SOINNs, the positive n-SOINN is considered to be the one with the low n value and the negative one with the high n value. The input vector of every SVM binary classifier is constructed by the output of the positive k n-SOINN and all other negative n-SOINNs. The output class along with its score of every binary SVM is then compared with all other binary SVMs in order to choose the top m classes. After choosing the top m classes, the output of their respective n-SOINNs is combined as an input to a multi-class SVM to get the final class. For our research, the validation was performed using the 41 attributes names of NSL-KDD dataset and all attributes information are available in [2]. The framework is based on the premise that TCP characteristics of a connection suggest whether a connection is defined as an attack or not and evermore if it is defined as an attack what kind of an attack that is.

2) Detection Engine.
The core module of the framework consists of two main parts as seen in Fig. 3: a clustering block and a classifier. The clustering block includes a pair of n-SOINNs there are used for each class to compress the information given from the TCP connections by the pre-processing module and achieve incremental learning. The classifying part takes the output of the n-SOINNs nodes constructs an input for an SVM classifier for each class to perform the preliminary classification. Subsequently, the top m classes of classification pairs of each classifier sorted by score are then classified by a multiclass SVM for the final decision (m is user defined according to the classes available). The score is based on the distance of the samples to the separating hyperplane.

3) Validation Module.
The validation module serves the purpose of improving and advance the accuracy of the framework by confirming the predicted label manually or automatically. The manual method could be a user interaction in real time or automatically with the predictions saved for later use and compared against known and confirmed predictions. After confirmation, the module forwards the failed predictions to the updating module.

4) Updating Module.
The framework operates in two phases: ,live and updating. In the live phase, it is making decisions based on what its capabilities where at that time and an update phase where the update module is updating the system with failed predictions to improve its capabilities. The phases could run in parallel if needed in production or alternating, according to the network traffic.

IV. EXPERIMENTAL RESULTS
Experimentation and evaluation were performed with the NSL-KDD dataset, which is an improved version of wellknown KDD'99 dataset, a very good candidate dataset to evaluate the performance of any IDS and thus our framework. For the negative SOINN we used n=2 and for the positive SOINN we used n=100. As the dataset included a mixture of normal traffic and four attacks, five binary classes were created for this experiment -a normal class and four attack classes (Denial of Service, Probe, R2L and U2R). For the top m classes, we choose m=3 (the smallest possible number higher than 2). The implementation used LIBSVM [24], which is the best performer in a multi-class SVM out of the box; we also used the RBF kernel in both binary and multiclass SVMs.
The dataset used included two different subsets, one with 125973 records and the other one with 22544 records. The second subset is not from the same probability distribution as the first one and it also contains specific types of attacks not present in the first one. For the framework evaluation and in order to show that incremental learning is achieved we have divided the first subset into five smaller subsets used for the test/update rounds and used the second subset with 22544 records for the initial training. It should be noted that after the initial training for every update round, the subset is tested against the trained n-SOINN-WTA-SVM and only the failed predictions are fed back to the system. The results, as shown in The accuracy chart, shown in Fig. 5, presents the relationship between the number of records accumulating and the accuracy of the system. The accuracy trend is pointing upwards with a 89.67% prediction result. The framework training samples accumulated by the end of Round 5 represented 27.30% of the total dataset; this makes it suitable for scaling applications because of its efficiency. Instead of feeding the complete dataset to the model, in our experimentation the SOINNs kept only a fraction of the information available for the whole dataset. In contrast, static models are trained with the full data available without the ability to adapt to the dynamic nature of network data over time. The time cost in seconds for the initial training was 986 seconds and after that, each update round the time cost incrementally increased with the last round costing 3285 seconds, as shown in Fig. 6. The fact that the offline method achieved less accurate results only make the case stronger that the incremental learning roposed framework could be a competitive candidate for an industrial application. The incremental learning property of the proposed framework and the evaluation results indicate that our proposal can adapt to the dynamic profile nature of network data for both normal and attacks categories. The framework utilizes less resources, is faster and it has higher detection rate than the offline method. By only feeding the system with failed predictions, not only we achieve incremental learning with promising accuracy but the framework is also efficient, meaning that, instead of feeding it the complete dataset, the SOINNs kept only a fraction of the dataset, making the framework a very good candidate for scaling systems.
Although the framework training time increases as the update data input grows, the framework update and live modes could either work in parallel (simultaneously) or the update mode could switch when the framework live phase is idle (i.e. no incoming data to detect).

VI. CONCLUSION AND FUTURE WORK
The security of the Internet of Things environments is the subject of research for almost two decades now and the rate of devices connected to the internet is growing exponentially. Each connected device is a potential target for an attack and IoT connected devices, excluding smartphones tablets and computers, have limited hardware resources which makes it even harder to apply proper security measures, increasing the exposure to vulnerabilities. The initial evaluation results shown that the proposed framework can achieve on-line updated incremental learning in a fast and efficient manner, making it suitable for an incrementally improve Network Intrusion Protection System. Our next steps are to improve the performance of the system, evaluate the proposed method with more recent datasets and attacks and finally to compare our results with state-of-the-art methods in terms of performance, False Positive and False Negative detection rate. Moreover, a possible combination of the technique with the emerging epidemic-based analysis to efficiently assess the probability of infection of a device in a network [25].