Johan Medrano
email: johan.medrano@lirmm.fr

Abderrahmane Kheddar
email: kheddar@lirmm.fr

Annick Lesne
email: annick.lesne@sorbonne-universite.frd

Sofiane Ramdani
email: sofiane.ramdani@umontpellier.fr

Radius selection using kernel density estimation for the computation of nonlinear measures

Keywords: Nonlinear measures, correlation sum, correlation dimension, Kolmogorov-Sinai entropy, kernel density estimation, recurrence plots

When nonlinear measures are estimated from sampled temporal signals with finitelength, a radius parameter must be carefully selected to avoid a poor estimation.

These measures are generally derived from the correlation integral which quantifies the probability of finding neighbors, i.e. pair of points spaced by less than the radius parameter. While each nonlinear measure comes with several specific empirical rules to select a radius value, we provide a systematic selection method.

We show that the optimal radius for nonlinear measures can be approximated by the optimal bandwidth of a Kernel Density Estimator (KDE) related to the correlation sum. The KDE framework provides non-parametric tools to approximate a density function from finite samples (e.g. histograms) and optimal methods to select a smoothing parameter, the bandwidth (e.g. bin width in histograms). We use results from KDE to derive a closed-form expression for the optimal radius. The latter is used to compute the correlation dimension and to construct recurrence plots yielding an estimate of Kolmogorov-Sinai entropy. We assess our method through numerical experiments on signals generated by nonlinear systems and experimental electroencephalographic time series.

Nonlinear measures computed from discrete time series are widely used to characterize and identify their underlying dynamics. Such measures include entropies, dimensions and indices that are derived from recurrence plots. These indices, all based on the fundamental concept of correlation sum, have been shown to be effective in distinguishing dynamical processes based on experimental data in various fields, e.g. mechanics, physiology, etc. One crucial parameter involved in their estimation is the radius used to define neighbors in the state space. Although many rule of thumbs are available, there is a need of well-grounded theoretical criteria to correctly select this radius parameter.

We propose to address this issue using Kernel Density Estimation (KDE). We first demonstrate the theoretical link between the correlation sum and the KDE framework. We then derive a loss function, whose optimization is equivalent to the minimization of the mean integrated squared error of a kernel density estimator, leading to a closed-form criterion for the radius parameter selection.

These findings moreover show how the estimator bias-variance trade-off determines a range for the radius values. Numerical experiments on both simulated (chaotic, corrupted by additive noise) and real-world data are presented to assess our approach.

I. INTRODUCTION

Nonlinearity and chaos govern a wide variety of systems. They are found in neurons firing patterns [START_REF] Faure | Is there chaos in the brain? I. Concepts of nonlinear dynamics and methods of investigation[END_REF]) and related electrophysiological signals [START_REF] Freeman | Evidence from human scalp electroencephalograms of global chaotic itinerancy[END_REF]), and in unpredictable changes of Earth climate [START_REF] Ghil | Climate dynamics and fluid mechanics: Natural variability and related uncertainties[END_REF]), to cite few examples.

Nonlinear measures of such systems are made more accurate thanks to an increasing interest in numerical tools suitable for nonlinear phenomena. Indeed, data generated by such systems are more suitable to nonlinear time series analysis, which provide complementary information to traditional linear methods such as power spectrum analysis [START_REF] Yang | A recurrence quantification analysis-based channel-frequency convolutional neural network for emotion recognition from EEG[END_REF]).

Our work focuses on estimating various metrics, measures (in the sense of quantitative indices), and invariants that rely on the computation of a correlation sum. The correlation sum is the estimator of the correlation integral, which is the mean probability that two points of the phase space trajectory of a dynamical system are neighbors [START_REF] Ott | Chaos in dynamical systems[END_REF]), i.e. the mean probability that their distance is less than a parameter called radius, threshold or tolerance depending on the application domain. The correlation sum captures important aspects of the nonlinear dynamics. Therefore, it is a fundamental quantity in various nonlinear measures: correlation dimension (Grassberger and Procaccia (1983a)), Kolmogorov-Sinai entropy [START_REF] Grassberger | Estimation of the Kolmogorov entropy from a chaotic signal[END_REF]; [START_REF] Eckmann | Ergodic theory of chaos and strange attractors[END_REF]; [START_REF] Faure | A new method to estimate the Kolmogorov entropy from recurrence plots: its application to neuronal signals[END_REF]), its approximate versions ApEn [START_REF] Pincus | Approximate entropy as a measure of system complexity[END_REF]) and SampEn [START_REF] Richman | Physiological time-series analysis using approximate entropy and sample entropy[END_REF]), Rényi's entropies [START_REF] Principe | Information theoretic learning: Renyi's entropy and kernel perspectives[END_REF]; [START_REF] Singh | Information theoretic learning with adaptive kernels[END_REF]), recurrence plots [START_REF] Eckmann | Recurrence plots of dynamical systems[END_REF], [START_REF] Marwan | Recurrence plots for the analysis of complex systems[END_REF]) and related metrics of recurrence quantification analysis [START_REF] Grendár | Strong laws for recurrence quantification analysis[END_REF]), etc.

In different nonlinear measures, the radius appears either as a variable or as a parameter. For instance, the correlation dimension is computed by estimating a scaling factor on a logarithmic plot of the correlation sum versus the radius. In contrast, a recurrence plot displays neighboring points on a black and white image and requires to fix the radius parameter beforehand. In both cases, the radius is selected as small as possible. As a correlation sum computed from a finite-length time series will likely tend to 0 together with the radius parameter, the challenge is to identify a radius range corresponding to a statistically useful distribution of neighbors. Eckmann and Ruelle [START_REF] Eckmann | Ergodic theory of chaos and strange attractors[END_REF]) (Section V.A.1.a.) refer to it as a "meaningful range" for the radius parameter. In our approach, we first derive an expression of the optimal radius. Then, we introduce a range to select a radius parameter or to study the properties of a function of the radius.

Several empirical rules exist to select a value or a range of values for the radius; however, they generally focused on a particular nonlinear measure [START_REF] Webber | Recurrence quantification analysis. Theory and Best Practices[END_REF]; [START_REF] Zbilut | Embeddings and delays as derived from quantification of recurrence plots[END_REF]; [START_REF] Pincus | Approximate entropy as a measure of system complexity[END_REF]). Here, we introduce a method which can be applied to any nonlinear measure derived from the correlation sum. Observing that logcorrelation sums are particularly used in nonlinear indices and that relative error arises from logarithmic error terms, we focus on minimizing a relative error between the correlation sum and the correlation integral. We show that minimizing the relative error term is equivalent to minimizing a well-known error used in the framework of Kernel Density Estimation (KDE), widely studied in statistics [START_REF] Silverman | Density estimation for statistics and data analysis[END_REF]) and in signal processing [START_REF] Gunduz | Correntropy as a novel measure for nonlinearity tests[END_REF]; [START_REF] Singh | Information theoretic learning with adaptive kernels[END_REF]). KDE denotes a family of non-parametric density estimation methods which generalize the well-known histogram methods [START_REF] Silverman | Density estimation for statistics and data analysis[END_REF]). Simple probability functions called kernels are placed at sample data points to approximate the underlying density function. In the KDE framework, the choice of the kernel width influences the degree of smoothing of the estimated density function. Selecting the kernel width is known as the bandwidth selection problem. The latter can be formulated simply as a bias-variance trade-off. Bandwidth selection is an extensively studied problem (see [START_REF] Jones | A brief survey of bandwidth selection for density estimation[END_REF] for a brief review) with notable usages in signal processing, e.g. mutual information estimation [START_REF] Moon | Estimation of mutual information using kernel density estimators[END_REF]).

The convergence of kernel density estimators for mixing dynamical systems was recently shown in [START_REF] Hang | Kernel density estimation for dynamical systems[END_REF].

The relation between kernel density estimation and the correlation sum is noted in [START_REF] Yu | Efficient implementation of the Gaussian kernel algorithm in estimating invariants and noise level from noisy time series data[END_REF] to estimate dynamical invariants in noisy situations. More recently, Gaussian kernels estimators of the correlation integral are applied to estimate Rényi's entropies [START_REF] Principe | Information theoretic learning: Renyi's entropy and kernel perspectives[END_REF]; [START_REF] Singh | Information theoretic learning with adaptive kernels[END_REF]; [START_REF] Erdogmus | From linear adaptive filtering to nonlinear information processing-the design and analysis of information processing systems[END_REF]). Here, KDE is used not to derive new estimators of nonlinear measures but rather as a framework providing a systematic rule to select the radius in computing nonlinear measures. We show that the radius minimizing the relative error of the correlation sum estimator is equivalent to the bandwidth minimizing the Mean Integrated Squared Error (MISE) of a density estimator (Section 3.1). Therefore, we use a bandwidth selection method from KDE to derive a closedform expression for the optimal radius (Section 3.2) and define a "meaningful range" for the radius variable relatively to our optimum (Section 3.3). We conduct numerical experiments on well-known dynamical systems. First, we study the behavior of the correlation sum estimator in the "meaningful range" for signals of different lengths and noise levels (Section 4).

Then, we estimate the Kolmogorov-Sinai entropy of both simulated and real signals, using recurrence plots computed with an optimal radius (Section 5).

II. CORRELATION SUM AND CORRELATION DIMENSION

Let (X , A, µ, T ) be a measure-preserving dynamical system with X ⊂ R d and µ the invariant measure (probability distribution in the phase space invariant upon the dynamics).

The correlation integral c(r) is the mean probability to find a pair of points at two different time x, y ∈ X arbitrarily close, such that the distance between x and y is less than a small radius parameter r [START_REF] Ott | Chaos in dynamical systems[END_REF]; [START_REF] Singh | Information theoretic learning with adaptive kernels[END_REF]):

c(r) = P ((x, y) : x -y p < r) = x∈X µ(B r (x))dµ(x) (1)
where B r (x) = {y ∈ X : x-y p < r} is the generalized d-dimensional ball in L p space, with radius r and center x. In practice, an estimator of the correlation integral can be computed from a sample trajectory Grassberger and Procaccia (1983a)):

x i ∈ R d , 1 ≤ i ≤ n (
C(r, n) = 1 n 2 n i,j=1 Θ(r -x i -x j p ) ---→ n→∞ c(r) (2)
where Θ is Heaviside step function and C(r, n) is called the correlation sum [START_REF] Pesin | Dimension theory in dynamical systems: contemporary views and applications[END_REF]; Grassberger and Procaccia (1983a)). For small values of r, the correlation integral grows as a power law:

c(r) ≈ const × r D 2 (3)
The quantity D 2 is called the correlation dimension.

III. KERNEL DENSITY ESTIMATION

A probability density function f may be estimated by placing smoothing kernels at each sample point. A smoothing kernel K is defined as a valid probability density function, which satisfies [START_REF] Silverman | Density estimation for statistics and data analysis[END_REF]):

K(u)du = 1 ∀u ∈ R, K(u) ≥ 0 (4)
Without loss of generality, we introduce a scaled version of the kernel with a L p norm and a scaling factor h > 0, K h (u) = h -d K(u/h), which is a valid kernel when K is a valid kernel.

A simple kernel is the uniform or boxcar kernel, which remains constant over a domain:

K h (u) = 1 τ p,d h d Θ(h -u p ) ( 5
)
where Θ is Heaviside step function and τ p,d is the volume of the unit ball defined by the norm p in a d-dimensional space (see Appendix A). Given samples

x i ∈ R d , 1 ≤ i ≤ n, distributed
according to a density f , a kernel density estimator of f is:

fh (x) = n -1 n i=1 K h (x -x i ) (6)
While kernel density estimators are consistent for i.i.d. samples, independence between consecutive samples cannot generally be assumed for dynamical systems. [START_REF] Hang | Kernel density estimation for dynamical systems[END_REF]) showed that kernel density estimators are also consistent for dynamical systems with mixing properties and weakly-continuous density function (more specifically C-mixing systems with pointwise α-Hölder controllable density, see [START_REF] Hang | Kernel density estimation for dynamical systems[END_REF] defs. 1 and 2).

The bandwidth parameter h determines the "width" of the kernels and consequently the degree of smoothing of the estimator. A plethora of methods exist to select the bandwidth parameter, see [START_REF] Jones | A brief survey of bandwidth selection for density estimation[END_REF]. Among existing bandwidth selection methods, minimizing the Asymptotic Mean Integrated Squared Error (AMISE), a Taylor expansion of the MISE

of the estimator E[ R d (f (x) -f (x)) 2 dx],
is appealing for practical applications as it allows to derive a closed-form expression of an approximately optimal bandwidth [START_REF] Silverman | Density estimation for statistics and data analysis[END_REF]):

h AMISE = W 1 (K) × d n × [W 2 (K)] 2 × W 1 (∇ 2 f ) 1/(d+4) (7)
where the functionals W i are defined as

W 1 (g) = R d g 2 (x)dx and W 2 (g) = R d x 2 1 g(x)
dx, where x 1 is a scalar component of x. Reference rules [START_REF] Silverman | Density estimation for statistics and data analysis[END_REF]; [START_REF] Scott | On optimal and data-based histograms[END_REF]) can be easily obtained by replacing the unknown quantity W 1 (∇ 2 f ) with the quantity computed using a reference distribution, generally a Gaussian distribution.

IV. A REFERENCE RULE FOR THE OPTIMAL RADIUS

To derive the expression of an optimal radius, we proceed as follows. First, we show that the radius minimizing the relative error of the correlation sum estimator is equivalent to the bandwidth minimizing the MISE of a particular density estimator. Second, we derive the closed-form expression of the radius minimizing the AMISE of the estimator. Finally, we identify a meaningful range to select a variable radius.

A. Criterion to select the radius

The correlation integral, c(r) = E µ [µ(B r (x))], is generally estimated by the correlation sum, Ĉ(r, n) = n -1 n i=1 μ(B r (x i )) (Eq. (1) and Eq. ( 2)). To obtain a good estimation of the correlation sum, we shall minimize the error between an estimator of the invariant measure of a ball , μ(B r (x)), and the true quantity µ(B r (x)). However, minimizing such error is not sufficient to provide a good estimation for small r: the scale of the error decreases with r and systematically leads to the trivial solution r = 0. Indeed, when r decreases, the absolute error decreases while the relative error is multiplied by a factor proportional to 1/r (Eq. ( 3)) and consequently blows up. Therefore, we want to find the radius minimizing a relative error criterion on μ(B r (x)). Let λ be the Lebesgue measure, such that λ(B r ) is the volume of a ball with radius r. We use the fact that µ(B r (•)) is proportional to r (see [START_REF] Eckmann | Ergodic theory of chaos and strange attractors[END_REF], Section V.A.) and consequently that µ(B r (•)) ∝ λ(B r ) to simplify the expression of the relative error and express the following local relative error:

L(r, x) = E µ(B r (x)) -μ(B r (x)) λ(B r ) 2 (8)
where the expectation is taken over samples used to construct the estimator. Given a fixed

r, µ(B r (•)) is a bounded function on X . We denote ρ the normalized density of µ(B r (•)), such that ρ(x) = µ (B r (x)) X µ(B r (x))dx = µ (B r (x)) λ(B r ) (9)
and, similarly, ρr (x) = μ(B r (x))/λ(B r ) the estimator of the normalized density ρ. After replacing in Eq. ( 8), we obtain:

L(r, x) = E (ρ(x) -ρr (x)) 2 (10)
Then, integrating Eq. ( 10) over possible values of x gives a global criterion to select the optimal radius r opt :

r opt = arg min r L(r) (11)
where

L(r) = X E (ρ(x) -ρr (x)) 2 dx (12)
With simple manipulations, we see that Eq. ( 12) is indeed the MISE between the estimator of the normalized density ρr and the true normalized density ρ [START_REF] Silverman | Density estimation for statistics and data analysis[END_REF]). Finally, ρr can be identified by replacing μ(B r (x)) with ρr in the expression of the correlation sum estimator (Eq. ( 2)):

C(r, n) = λ(B r ) n n i=1 ρr (x i ) = 1 n 2 n i,j=1 Θ(r -x i -x j ) yielding ρr (x) = 1 nλ(B r ) n i=1 Θ(r -x -x i ) (13)
As λ(B p,d r ) = λ(B p,d 1 )r d (A), we observe that ρr is a kernel density estimator with a uniform kernel and a bandwidth parameter r (Eq. ( 5)). Hence, it follows from Eq. ( 11) that the bandwidth minimizing the MISE of the estimator ρr can provide a good approximation of the optimal radius r opt minimizing the relative error on the correlation sum estimator.

Moreover, the AMISE method (Eq. ( 7)) can be used to approximate r opt with a simple, closed-form expression that resembles to the empirical rules currently used. In the next section, we use the AMISE minimization method to derive a reference rule for the optimal radius.

B. Derivation of a reference rule for the optimal radius

As presented in section III, a Taylor expansion of the MISE can be used to derive a closedform expression of the optimal bandwidth for an estimator. A particular interest for this method is motivated by the possibility of deriving a closed-form expression of the optimal bandwidth. We use a reference Gaussian distribution in Eq. ( 7) and derive the expressions for W 1 (K) and W 2 (K) for the uniform kernel (see Appendix B). Then, substituting these expressions into Eq. ( 7) gives the main result of the paper: a reference rule radius r opt defined as

r opt = α p,d × ŝ × n -1/(d+4) (14)
where α p,d , depending on the norm and dimension, rescales r opt ; ŝ is an estimate of the spread of data; and n is the length of the trajectory in phase space.

Remark 1. In practice, the phase space is reconstructed using a time delay embedding procedure (according to Takens theorem [START_REF] Takens | Detecting strange attractors in turbulence[END_REF])) ; hence, if N denotes the length of the univariate time series, d the embedding dimension and τ the delay, the length of the trajectory in reconstructed phase space is n = N -(d -1)τ .

Estimation of the spread ŝ

A first choice for the spread ŝ is the average marginal sample standard deviation, defined by σ = d -1 i S i,i , S ∈ R d×d is the sample covariance matrix. When the d-dimensional sample is constructed from an univariate time series using delay embedding, components on the diagonal of the sample covariance matrix are equal: σ is then the sample standard deviation of the time series. Alternatively, as the interquartile range IQR is a good alternative to standard deviation for non-Gaussian data (see [START_REF] Silverman | Density estimation for statistics and data analysis[END_REF] for discussion), a common choice for ŝ is:

ŝ = min σ, IQR 1.34 . ( 15
)
2. Derivation of the reference factor α p,d

The expression for the 1-dimensional reference factor is relatively straightforward:

α p,1 = (12 √ π) 1/5 ≈ 1.843.
The general closed-form expression for α p,d is more complex (see Appendix B, Eq. ( B7)); however, the expression can be simplified for common norms (Appendix C):

α 1,d = (d + 2)! (d + 1)( √ π) d 1/(d+4) (16) α 2,d = 2 × Γ d 2 + 2 2 1/(d+4) (17) α ∞,d = 36( √ π) d d + 2 1/(d+4) (18)
Moreover, α p,d is to be computed only once for common dimensions and norms. Hence, we report in Table I some values of α p,d that can be used in Eq. ( 14).

C. Identification of a meaningful range for a variable radius

As discussed above, some nonlinear indices require selecting the range of radius values in which the quantity is estimated. For instance, this applies to nonlinear indices quantifying a scaling exponent of the form lim r→0 log ν(Br) log r

(with ν a (with nu a probability distribution in the phase space), as often encountered in the chaotic systems literature (see e.g. [START_REF] Pesin | Dimension theory in dynamical systems: contemporary views and applications[END_REF];

Ott ( 2002)). In practice, the limit r → 0 is generally intractable, and estimations of ν for

α p,d p 1 2 ∞ d 1 1.843 1.
r opt = α p,d × min (σ, IQR/1.34) × n -1/(d+4) .
small r are highly variable due to poor statistics [START_REF] Eckmann | Ergodic theory of chaos and strange attractors[END_REF]). On the other hand, at a certain point, large r will not capture the desired scaling effect. Hence, there is a range of values which must be selected to support a good estimation of the nonlinear measure. Here, we introduce our arguments to guide the selection of a meaningful range for a variable radius.

The AMISE can be expanded in an integrated squared bias and integrated variance of the density estimator, giving the following expressions of bias and variance as functions of r (Silverman (1986)):

bias(r) r 2 2 W 2 (K)∇ 2 ρ(x) (19) var(r) W 1 (K) nr d ρ(x) (20)
The behavior of the relative error with r can be understood from Eq. ( 19) and Eq. ( 20):

the bias is proportional to r whereas the variance is inversely proportional to r. As r opt minimizes the AMISE, the bias contribution increases with r while the variance decreases with r. However, the bias term only depends on r whereas the variance term decreases when the number of points increases. This observation -considering that usually the bandwidth minimizing the AMISE is too large, suggests selecting the optimal radius r opt as the upper bound for the meaningful range. We introduce a range parameter 0 < β < 1 to select the lower bound as a fraction of r opt , such that the radius values lie within the range:

R = [βr opt , r opt ] (21)
Due to the relations Eq. ( 19) and Eq. ( 20), we argue that the value of β shall be decreased when increasing the number of points.

V. ESTIMATION OF THE CORRELATION DIMENSION

In the following, we investigate the behavior of the Grassberger and Proccacia algorithm for the estimation of the correlation dimension. We compare the spread and bias of estimations in the full range of available scales with estimations in the meaningful range derived in section IV C.

A. The Grassberger and Proccacia algorithm

The correlation dimension D 2 can be expressed as:

D 2 = lim r→0 log c(r) log r = lim r→0 lim n→∞ log C(r, n) log r (22)
The Grassberger and Proccacia algorithm (Grassberger and Procaccia (1983a)) for the empirical estimation of the correlation dimension consists in computing the correlation sum for different values of r and plotting log C(r, n) versus log r. The slope of the linear region in this logarithmic plot provide the desired estimation of the correlation dimension D 2 [START_REF] Ott | Chaos in dynamical systems[END_REF]). Similarly to the original paper (Grassberger and Procaccia (1983a)), we use linear regression to estimate the slope.

B. Procedure for generating reconstructed trajectories

We conducted numerical experiments on the Lorenz system [START_REF] Lorenz | Deterministic nonperiodic flow[END_REF]) (σ = 10, β = 8 3 , ρ = 28, dt = 0.01), the Rössler system [START_REF] Rössler | An equation for continuous chaos[END_REF]) (a = 0.1, b = 0.1, c = 14, dt = 0.05) and the Hénon map [START_REF] Hénon | A two-dimensional mapping with a strange attractor[END_REF]) (a = 1.4, b = 0.3). We apply the following procedure to generate random time series of different length. After drawing a random initial state, we generate time series for all systems -using a Runge-Kutta 4/5 method for Lorenz and Rössler -such that the length of the time series is N after removing transients (sample series are presented in Figure 1). Then, the trajectory is reconstructed using Takens delay embedding (series from the x coordinates were systematically used). The original system dimension is used as embedding dimension d. The time delay parameter τ is set to 1 for the Hénon map and selected as the first minimum of the time-delayed mutual information function [START_REF] Fraser | Independent coordinates for strange attractors from mutual information[END_REF]) for the Lorenz and Rössler systems. Please note that the length of the reconstructed trajectories, n = N -(d-1)τ , is used to compute the optimal radius using Eq. ( 14) (see Remark 1).

Remark 2. Here, we assume that the delay and embedding dimension are correctly selected as a bad phase space reconstruction deteriorates the log C(r, n) versus log r plot [START_REF] Kantz | Nonlinear time series analysis[END_REF]). In practice, this embedding problem can be efficiently addressed as a plethora of methods exist to select the delay and the embedding dimension (see for instance [START_REF] Kantz | Nonlinear time series analysis[END_REF], Ch. 3.3 and Ch. 9.2).

C. Numerical results for the radius range

Here, we visualize the meaningful range on the log-log plot of the correlation sum versus the radius. We first generate 100 series for each system (4000 points for Rössler and Lorenz attractor, 200 points for Hénon map, respectively), select 25 random values of radius and compute the corresponding correlation sums. We overlay the average value of r opt and the ranges with arbitrary values β ∈ {0.01, 0.1, 0.5} on the plot of log C(r, n) vs log r. Results are presented in Figure 2. We observe that the spread of the correlation sum over the runs is low at the location of reference radius and increases when the radius is decreased. Hence, smaller values of β likely lead to higher variance estimations.

On Figure 2b, a knee is present around a value log(r knee ) -1, such that the slope a left in a left range [r 0 , r knee ] is higher than the slope a right in right range [r knee , r 1 ]. In practice, a knee may appear from the superposition of signals from non-interacting subsystems with different amplitude [START_REF] Eckmann | Ergodic theory of chaos and strange attractors[END_REF]). In this situation, a left characterizes the two subsystems while a right corresponds only to the system with the largest signal amplitude.

Consequently, a careful analysis of the plot of log C(r, n) vs log r might be necessary to select a range capturing the desired properties of systems under study.

D. Influence of the time series length

Using the procedure described in Section V B, we generate 100 trajectories for each length:

(a) N = 250, 500, 1000, 2500, 5000 for the Lorenz system, (b) N = 500, 1000, 2500, 5000, 7500 for the Rössler system, (c) N = 100, 250, 500, 1000, 2500 for the Hénon map.

Remark 3. Notice that the discrepancies for the number of points used for the three systems can be justified by the resulting trajectories after time delay embedding. Indeed, when the number of points is too low, the reconstructed trajectories cannot properly reflect the dynamics nor the correct dimension of the attractor. For instance, in our experiments this was the case for time series of 100 points for the Rössler system.

We computed correlations sums for 20 values of r ranging between 10 -8 and 2σ, where σ is the sample standard deviation. We compared the estimation using the Grassberger and Proccacia algorithm on the entire curve (the plateau on the right was omitted) with the estimation suggests that the range must be selected sufficiently large to provide a proper support for dimension estimation. Moreover, although the bias of the Grassberger-Proccacia algorithm is low in this setup, a single dimension estimate can be far from the true dimension. Therefore, one can favor a smaller range for r to reliably estimate a quantity slightly lower than the true dimension.

We found qualitatively similar results for Lorenz and Rössler attractors when series of different length are obtained by downsampling an original series of fixed length (results not shown).

E. Influence of observational white noise

Finally, we investigate the influence of observational noise on the estimation of the correlation dimension in the different ranges. Observational noise is ubiquitous in practical applications and creates a knee on the plot of log C(r, n) versus log r, with a dimension at the left of the knee equal to the embedding dimension (see [START_REF] Eckmann | Ergodic theory of chaos and strange attractors[END_REF]; [START_REF] Grassberger | Measuring the strangeness of strange attractors[END_REF]). Hence, the range must be selected at the right of the knee to provide good estimations of the dimension.

We generate 100 time series of 1000 points for the three systems. Each series, with standard deviation σ, is corrupted with additive white Gaussian noise with standard deviation σ noise = k σ, where k defines the noise level. As above, we compare the estimation of the original Grassberger and Proccacia algorithm with the estimation in the different ranges (the reference radius Eq. ( 14) is computed for each noise-corrupted series). We present in original version of the Grassberger and Proccacia algorithm (red) or the range β = 0.01 (beige). In contrast, larger values of β yield more consistent results under the different noise conditions. Therefore, this observation suggests that in noise conditions, the correlation dimension can be more robustly estimated from a smaller range of r.

VI. ESTIMATION OF KOLMOGOROV-SINAI ENTROPY USING RECURRENCE PLOTS

In this section, we investigate the behavior of the goodness-of-fit of an estimator of a nonlinear measure with the radius parameter. We also study the reference radius inherent to Eq. ( 14) under different conditions. We use the reference radius in the construction of recurrence plots used to estimate the Kolmogorov-Sinai (KS) entropy of Hénon map and apply similar method to real electroencephalographic (EEG) signals.

A. Recurrence plots and Kolmogorov-Sinai entropy

Recurrence plots

Recurrence plots [START_REF] Eckmann | Recurrence plots of dynamical systems[END_REF]) display phase-space neighbors as a 2D blackand-white image whose (i, j) element is black if trajectory points x i and x j are closer than a fixed radius ε. More formally, from a phase-space trajectory {x i }, 1 ≤ i ≤ n, a recurrence plot rp(ε) ∈ R n×n is defined as:

(rp(ε)) i,j = Θ(ε -x i -x j p ) ( 23
)
where

Θ(•) denotes Heaviside step function, • p is a norm, usually either L 1 , L 2 , or L ∞ .
The patterns in recurrence plots reflect properties of the underlying dynamical system and can be quantified using the Recurrence Quantification Analysis (RQA) framework, providing a set of powerful non-parametric visualization and characterization tools for nonlinear time series analysis. The relationship between recurrence plots (and RQA measures) and the correlation sum intuitively follows Eq. ( 23) [START_REF] Grendár | Strong laws for recurrence quantification analysis[END_REF]). Indeed, simple mathematical manipulations show that the recurrence rate, defined as the average number of recurrent points in a recurrence plot, is equal to the correlation sum [START_REF] Thiel | Analytical description of recurrence plots of white noise and chaotic processes[END_REF]).

Estimating Kolmogorov-Sinai entropy from recurrence plots

The Kolmogorov-Sinai (KS) or measure-theoretic entropy [START_REF] Kolmogorov | A new metric invariant of transitive dynamical systems and automorphisms of Lebesgue spaces[END_REF], Sinai Ya (1959)) measures the evolution of uncertainty with the iteration of the map of a dynamical system. The lower bound K 2 , often used as the estimate of KS entropy [START_REF] Faure | A new method to estimate the Kolmogorov entropy from recurrence plots: its application to neuronal signals[END_REF]), is defined as [START_REF] Grassberger | Estimation of the Kolmogorov entropy from a chaotic signal[END_REF]):

K 2 = lim r→0 lim m→∞ lim n→∞ 1 ∆t log C m (r, n) C m+1 (r, n) (24)
C m (r, n) denotes the correlation sum built from a delay-reconstructed trajectory in a mdimensional L ∞ space. While it is possible to approximate the KS entropy directly from correlation sums [START_REF] Pincus | Approximate entropy as a measure of system complexity[END_REF]; [START_REF] Richman | Physiological time-series analysis using approximate entropy and sample entropy[END_REF]), we rather consider the method in [START_REF] Faure | A new method to estimate the Kolmogorov entropy from recurrence plots: its application to neuronal signals[END_REF]). The latter approximates the KS entropy from the histogram of diagonal lines of length greater than m in a recurrence plot rp(ε):

N ε (m) = card{(i, j) : ∀k ∈ {0, . . . , m -1}, |u i+k -u j+k | < ε} (25)
A diagonal of size m on a recurrence plot reflects that two trajectories stayed at a distance smaller than a threshold ε for m time-steps, or equivalently that two delay-reconstructed vectors in m-dimensional space are close under L ∞ norm. Hence, the histogram of diagonal lines, N ε (m), captures information similar to the correlation sum from delay-coordinates, C m (r, n); whereas the parameters ε and r are analogous in the two quantities. The main advantage of the Faure and Korn method is computational: while C m (r, n) is computed for several values of the embedding dimension m, the histogram N ε (m) is computed only once.

Then, using the diagonal line histograms to rewrite the KS entropy Eq. ( 24) as a function of r = ε gives:

K 2 (r) = lim m→∞ lim n→∞ 1 ∆t log N r (m + 1) N r (m + 2) (26)
Faure and Korn [START_REF] Faure | A new method to estimate the Kolmogorov entropy from recurrence plots: its application to neuronal signals[END_REF]) suggest to evaluate the average slope of a log N r (m) vs m plot for various values of r. Then, taking the limit r → 0 is supposed to converge to a constant value equal to the KS entropy, K 2 , up to a scaling factor. However, selecting the smallest possible r to estimate the limit r → 0 from real-world data (i.e. finite-size samples with noise) likely leads to estimations flawed by a large variance, as discussed in [START_REF] Faure | Estimating Kolmogorov entropy from recurrence plots[END_REF]. Hence, the problem is to select a value of r yielding the best possible estimations of the KS entropy. We use our reference radius (Eq. ( 14)) to compute the recurrence plot used to estimate the KS entropy. Recurrence plots and diagonal line histograms were computed using the pyunicorn package [START_REF] Donges | Unified functional network and nonlinear time series analysis for complex systems science: The pyunicorn package[END_REF]).

B. Numerical experiments for the Hénon map

We generate 100 time series for each length (n = 150, 250, 500, 1500 points) from the standard Hénon map (a = 1.4, b = 0.3). For each series, we use the Faure and Korn method and compute the value of K 2 (r), Eq. ( 26), as a function log r curve for 50 values of log r ranging from -4 to 0.5. We then compute the reference radius (Eq. ( 14)) -using the series length n and dimension d = 1 -and average the values over series of same length. Results are presented in Figure 5. We notice that the variance of the estimation increases for decreasing radius and decreasing number of points. This result is presumably due to a poor statistical power for small values of the radius and short time series. However, the estimation seems to converge in average to the theoretical value (H KS = 0.42 (Faure and Korn (1998))) when r tends to 0.

Notice that the right-most part of the plot exhibits a large bias between the estimated and theoretical entropy values. Contrary to the variance, the bias does not seem to decrease with increasing number of points. Thus this bias is more symptomatic of the radius being too high to obtain any valuable information about the Hénon map. This bias-variance trade-off is usually related to a Mean Squared Error (MSE) minimization problem.

-

4 -3 -2 -1 0 -7 -6 -5 -4 -3 -2 n = 150 n = 250 n = 500 n = 1500
FIG. 6: Estimation of the log-Mean Square Error of the Kolmogorov-Sinai entropy estimator as a function of the radius value (lower is better). The filled areas corresponds to the 95% bootstrap confidence interval for each length. We see that for different time series, the reference radius gives a log-MSE value between -7 and -5.

As the MSE of an estimator quantifies the goodness-of-fit, the parameters of the estimator yielding the minimum value of MSE can be systematically selected. We numerically compute the MSE of the KS entropy estimator as a function of log r and use this plot as an objective criterion to evaluate the adequacy of our reference radius. The MSE consists in the sum of a squared bias term, measuring the difference between the theoretical value and the estimation, as well as the variance of the estimator. We use a theoretical value K 2 = 0.42 and all of the 100 sample series to compute the MSE, overlay the reference radius averaged over series of same length, and show the results in Figure 6. Each record is divided in four segments of 1024 points. For each segment, we compute a recurrence plot with the radius set by Eq. ( 14) and estimate the KS entropy using the for the healthy volunteers (control group) and the epileptic patients. Our estimator gives an average KS entropy of 0.288 ± 0.005 (95% confidence interval) for the epileptic group and 0.504 ± 0.006 for the control group, which confirms an average significant decrease of the KS entropy with epilepsy, as reported in previous studies [START_REF] Kannathal | Entropies for detection of epilepsy in eeg[END_REF]).

Finally, to compare the discrimination strength of common closed-form radius selection methods, we estimate the KS entropy with each method, perform a two-samples Z-test (epileptic versus control group) and collect the Z-score. We report a Z-score of Z = 45.3 (resp. Z = 39.3) for the r = 0.2σ (resp. r = 0.1σ, with σ the series standard deviation) rule [START_REF] Pincus | Approximate entropy as a measure of system complexity[END_REF]), Z = 48.9 when the radius is set to 10% of the maximum phase space neighborhood [START_REF] Zbilut | Embeddings and delays as derived from quantification of recurrence plots[END_REF]), Z = 41.1 (resp. Z = 34.6) when the radius is selected such that 10% (resp. 4%) of the number of points are selected as neighbors plots for an healthy volunteer and an epileptic patient respectively. The value r = 1.843 × ŝ × n -1/5 (Eq. ( 14)) is used to compute the recurrence plots from the univariate time series; (c) contains a box plot of the estimated entropies for both control and epileptic groups. As expected, the entropy values are significantly lower for the epileptic group. [START_REF] Webber | Recurrence quantification analysis. Theory and Best Practices[END_REF], [START_REF] Kraemer | Recurrence threshold selection for obtaining robust recurrence characteristics in different embedding dimensions[END_REF]), Z = 54.9 for the reference rule radius r = 1.843 × ŝ × n -1/5 (Eq. ( 14)). Subsequently, although all methods detect significant differences between the two groups, the radius given by Eq. ( 14) gives the most statistically significant results.

VII. DISCUSSION AND CONCLUSION

We propose a new approach for selecting the radius parameter in nonlinear measures derived from the correlation sum. We first formulate a relative error function on the quantities underlying correlation sums. We show that minimizing the loss function is equivalent to minimizing the MISE of a kernel density estimator. We use the AMISE minimization method to derive a closed-form expression to select the radius. Additionally, we observe how the bias and variance of the estimator varies with the radius and derive a "meaningful" range to select a variable radius.

We investigate the behavior of the Grassberger and Proccacia algorithm for estimating the correlation dimension in radius ranges of different size. We observe that the range parameter β can be selected close to 1 for low-variance estimations, and close to 0 for lowbias estimations. However, the presence of noise in the observed signal induces typical error in the estimations and leads to favor small ranges close to the reference radius.

We then use the reference radius to construct recurrence plots for estimating the Kolmogorov-Sinai entropy from both simulated and experimental signals. In a first analysis, we reconstruct the Mean Squared Error curve of the entropy estimator for Hénon map and show that the reference radius is close to the minimum of the curve. We confirm the experimental adequacy of the method by obtaining significant results in characterizing epileptic EEG signals.

Moreover, our theoretical approach yields a reference radius that is similar to several existing radius selection methods arising from empirical or numerical experience: the radius is a fraction of the scale of the data [START_REF] Pincus | Approximate entropy as a measure of system complexity[END_REF], [START_REF] Zbilut | Embeddings and delays as derived from quantification of recurrence plots[END_REF]) and compensates for the dimension of the data [START_REF] Kraemer | Recurrence threshold selection for obtaining robust recurrence characteristics in different embedding dimensions[END_REF]).

For the specific case of recurrence plots, [START_REF] Andreadis | On a topological criterion to select a recurrence threshold[END_REF] recently proposed an empirical procedure to identify an optimal radius value. They define a metric to measure the distance between recurrence plots and compute the distance between recurrence plots constructed from the same time series using increasing values of radius. The radius value is considered "optimal" when it minimizes the distance between consecutive recurrence plots, i.e. such that a slightly changing the radius has the minimal impact on the recurrence plot.

The principal issues with this procedure are the computational burden of building several recurrence plots and the difficulty to reliably identify the optimum. In contrast, our method is computationally much more efficient and not restricted to recurrence plots.

Our numerical experiments suggest that the reference radius given in Eq. ( 14) can be used as a default parameter to obtain robust and significant values for a number of different nonlinear tools and measures: correlation dimension, recurrence plots, Kolmogorov-Sinai entropy. In future work, we plan to investigate the relation between our optimal radius and the embedding parameters, which play a role on the trajectories resolution in the reconstructed phase space. Additionally, we plan to use the reference radius in EEG signal processing application, notably to extract dynamical features characterizing the oscillatory The general formula for d ≥ 2 is more complex:

r opt =    4d(2 √ π) 2 τ p,d d(d + 2) × 3Γ d+2 p + 1 Γ 1 p + 1 2 Γ d p + 1 Γ 1 + 3 p 2    1/(d+4) • ŝ • n -1/(d+4) =    4(2 √ π) d 3Γ d+2 p + 1 Γ 1 p + 1 2 τ p,d (d + 2) Γ d p + 1 Γ 1 + 3 p 2    1/(d+4) • ŝ • n -1/(d+4) (

FIG. 1 :

 1 FIG. 1: Sample series for the Lorenz attractor (a), the Rössler attractor (b) and the Hénon map (c). Series from the x coordinates were systematically used to reconstruct the trajectory using Takens delay embedding.

  FIG. 2: Log-log plot of correlation sums versus the radius used to estimate the correlation dimension for Lorenz (a), Rössler (b) and Hénon (c) systems. The dots correspond to estimates of the correlation sum at random values of radius for trajectories integrated from random initial states. The colored regions show different ranges defined by [βr opt , r opt ], with β = 0.01 (in beige), 0.1 (in light blue) or 0.5 (in dark blue).

  FIG. 3: Influence of the number of points on the estimation of the correlation dimension from Lorenz (a), Rössler (b) and Hénon (c) systems. We compare the original version of the Grassberger and Proccacia algorithm (red) with an estimation of the slope in the range [βr opt , r opt ], with β = 0.01 (beige), β = 0.1 (light blue), β = 0.5 (dark blue).

Figure 4 a

 4 Figure 4 a violin plot for noise levels k = 0, 0.05, 0.1, 0.15, 0.2. For both Rössler and Lorenz attractors, we observe that a noise level of 5% is sufficient to corrupt estimations with the

FIG. 5 :

 5 FIG. 5: Estimation of Kolmogorov-Sinai entropy for time series from the Hénon map with different lengths n. The filled areas corresponds to the 95% (Gaussian) confidence intervals for each length. The vertical dashed lines represents the average reference radius associated to each length. The horizontal dashed line indicates the reported entropy for Hénon map, K 2 = 0.42.

  For short time series, we observe that the radius selected by the reference rule is systematically close to the minimum of the MSE. For longer time series, the reference radius is larger than the minimum of the curve. Nevertheless, for values of r < r opt , the slope of the MSE curve gets flatter for increasing number of points and allows arbitrary selection of smaller radius values. We report similar observations for two other estimators of the KS entropy, the Approximate and Sample entropies (results not shown).C. Application to EEG signals in the context of epilepsyTo show the viability of our approach on real-world data, we apply our radius selection procedure to estimate the KS entropy of epileptic EEG signals. A significant decrease of the EEG signal entropy at the epileptic seizure location is a common feature for automatic seizure detection[START_REF] Ocak | Optimal classification of epileptic seizures in EEG using wavelet analysis and genetic algorithm[END_REF];[START_REF] Srinivasan | Approximate entropy-based epileptic EEG detection using artificial neural networks[END_REF]). We use the data publicly available from the University of Bonn[START_REF] Andrzejak | Indications of nonlinear deterministic and finite-dimensional structures in time series of brain electrical activity: Dependence on recording region and brain state[END_REF]), which consists in five sets of EEG data. Each set contains 100 segments of 23.6 seconds recorded at 173.61Hz (4096 points per segment), which were visually inspected for artifacts and band-pass filtered between 0.5Hz and 40Hz. Two sets contains surface EEG recorded from five healthy volunteers at rest, either with closed (set O) and opened eyes (set Z). The three other sets, consisting in signals from five epileptic patients recorded during presurgical evaluation, contain segments either from seizure-free intervals (at epileptogenic site, set F, or at the hippocampal formation of the opposite hemisphere of the brain, set N) or during seizure (at epileptogenic site, set S).

Faure

  and Korn method. Recurrence plots and signals sampled from the sets Z and F are shown in Figure 7a and Figure 7b. We present in Figure 7c a box plot of the KS entropy

  FIG. 7: Estimation of the Kolmogorov-Sinai entropy from recurrence plots to discriminate epileptic from healthy EEG signals: (a) and (b) show sample EEG signals and recurrence

  This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. Grassberger-Proccacia β = 0.01 β = 0.1 β = 0.5Noise level (%) Correlation dimensionThis is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0055797

TABLE I :

 I Rounded values of the coefficient α p,d for common norms and dimensions. The values can directly be used in reference rule radius,

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0055797

ACKNOWLEDGMENTS

This research was supported by a grant from the French Ministère de l'Éducation Nationale, de l'Enseignement Supérieur et de la Recherche.

DATA AVAILABILITY

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

dynamics of motor imagery EEG signals.

Appendix A: Volume of generalized balls Let B p,d r (x) be an open d-ball of size r in an L p space, i.e. B p,d r (x) = {y ∈ R d : x -y p < r}. We write τ p,d = λ(B p,d 1 ) the volume of the generalized unit ball, where λ is the Lebesgue measure. The volume of a generalized ball of radius r is λ(B p,d r ) = τ p,d r d . The general formula for the volume of a generalized unit ball is [START_REF] Wang | Volumes of generalized unit balls[END_REF]):

which can be simplified for common L p spaces:

Γ denotes Euler's gamma function with the property Γ(z + 1) = zΓ(z).

Appendix B: Derivation of a reference rule for the uniform kernel

The expression of the bandwidth minimizing the Asymptotic Mean Integrated Squared

Error is [START_REF] Silverman | Density estimation for statistics and data analysis[END_REF], Eq. 4.14 and 4.15):

where W i are the functionals W 1 (g) = R d g 2 (u)du and W 2 (g) = R d u 2 1 g(u)du, where u 1 is the first component of u ∈ R d (as the kernel is symmetric, it is sufficient to consider only u 1 in W 2 ). We compute W 1 for the uniform kernel:

W 1 (∇ 2 ρ) for a d-dimensional Gaussian reference distribution φ is given in (Silverman, 1986, Eq. 4.13):

Then, W 2 in the 1-dimensional case:

Changing to spherical coordinates u 2:d = (η, ξ) with an orientation vector ξ and a radius η = ( d i=2 |u i | p ) 1/p and :

where 2 F 1 (a, b; c; z) is the Gaussian hypergeometric function. Finally, using the fact that

) , (see Appendix A.1) and the expansion of 2 F 1 at z = 1 (Olver et al., 2020, Eq. 15.4.20), we can further simplify:

We then derive the reference rule by plugging the appropriate values in Eq. (B1). First, we derive the expression in the simple 1 dimensional case, which is independent from p:

Appendix C: Simplification of the reference rule for common norms

We address the case p = 1:

Then, the limiting case p → ∞:

Finally the case p = 2:

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS