# On Consistency of Cross-Approximate Entropy in Cardiovascular and Artificial Environments

Authors: Tamara Skoric; Olivera Sarenac; Branislav Milovanovic; Nina Japundzic-Zigon; Dragana Bajic
Journal: Complexity
 (2017)

---

## Abstract

Cross-approximate entropy(XApEn) quantifies the mutual orderliness of simultaneously recorded time series. Despite being derived from the firmly established solitary entropies, it has never reached their reputation and deployment. The aim of this study is to identify the problems that preclude widerXApEn implementation and to develop a set of solutions. Exact expressions forXApEn and its constitutive parts are derived and compared to values estimated from artificial data. This comparison revealed vast regions within the parameter space that do not guarantee reliable probability estimation, makingXApEn estimates inconsistent. A simple correction to one of theXApEn procedural steps is proposed. Three sets of formulae for joint parameter selection are derived. The first one is intended to maximize threshold profile. The remaining ones minimizeXApEn instability according to the strong and weak criteria. The derived expressions are verified using cardiovascular signals recorded from rats (long signals) and healthy volunteers (short clinical signals), proposing a change of traditional parameter guidelines.


---

## Body

## 1. Introduction

Approximate entropy(ApEn) is among the most exploited nonlinear techniques to quantify the complexity and unpredictability of time series [1–3], with a citation rate highlighted in [4].ApEn is approved as a supporting tool in preclinical and clinical studies, with the most prominent applications in cardiovascular studies.In spite of the vast and firmly established implementation, the activities onApEn improvement have never ceased. Research efforts are oriented towards modifications, including sample entropy(SampEn) [5], corrected conditional entropy(CApEn) [6], K-nearest neighborhood entropy(KNNCE) [7], fuzzy entropy(FuzzyEn) [8], multiscale entropy(MultiScaleEn) [9], distance entropy(DistEn) [10], binarized entropy(BinEn) [11],ApEn based on wave mode [12], and speeding-up algorithms [13, 14], and towards the consistency studies, including threshold selection [15–20], time lag [21], and consistency in general [4].One of the very firstApEn modifications isCross-ApEn (XApEn) [22] that estimates mutual predictability of two simultaneously recorded time series. In view ofApEn as a complexity measure of solitary signal,XApEn should have become an equivalent tool for paired physiological processes. However, the reputation and deployment of its predecessor have never been reached. Except for a recent improvement [23],XApEn applications follow classicalApEn guidelines, although (a) the parameters should be adjusted even forApEn (e.g., [4, 15–20]) and (b) properties ofXApEn are different and straightforwardApEn guidelines might not be proper.The aim of this paper is to identify the problems that preclude widerXApEn implementation and to propose methods for their adjustment. It is proved that the fundamental source ofXApEn inconsistency is unreliable estimation of conditional probabilities that are its constitutive elements. It is shown that the parameters required forXApEn are mutually dependent, so the traditional guidelines for their choice need an update. In order to locate the unstable parameter regions, theXApEn values estimated from the artificially simulated data are compared to the trueXApEn values calculated by a formula. Following these observations, we derived equations that jointly select stable parameters for reliableXApEn estimation in real signals, for which exact mathematical descriptions are unavailable. The derived equations are tested using a wide variety of cardiovascular signals: using typical short records of ambulatory patients and using long signals recorded from small laboratory animals.


## 2. Materials and Methods

### 2.1. Experimental Data

The study used three types of data: artificial, animal, and human. The purpose of artificial data was to test the stability of estimated values, comparing them to the exact values calculated by the derived formulae, and also to develop a new set of threshold equations. The data recorded from animals were also used for method development and for its repeated testing, enabled by the record length. The data recorded from human volunteers were used to validate the derived methods in typical real applications where the records are short.The animal data include systolic blood pressure (SBP) signals and pulse interval (PI) signals recorded from two groups of rats: the first group comprised outbred male Wistar normotensive rats, while the second group comprised outbred male borderline hypertensive rats. Ten days before the recordings, radiotelemetric probes (TA11PA-C40, DSI, Transoma Medical) were implanted in the abdominal aorta. Blood pressure (BP) waveforms were digitally recorded and transmitted using Dataquest A.R.T. 4.0 equipment (DSI, St. Paul, MN, USA) with sampling frequencyfs = 1000 Hz. The systolic blood pressure (SBP) was determined at the receiver as waveform local maxima. The pulse interval was determined as an interval between the successive moments of maximal gradients ΔBP/ΔTmax, that is, the local maximal change (Δ) of blood pressure (BP) waveform per sampling interval ΔT. Pulse interval is an acceptable alternative to R-R interval of an electrocardiogram. During the experiments, the animals were exposed to two types of stress: shaker stress, with rats positioned on a platform shaking at 200 cycles/min, and restraint stress, with rats placed in a Plexiglas restrainer tube (ID 5.5 cm with pores) in the supine position, with a detailed description given in [24]. These data were already used for other research purposes [20, 24, 25]. The duration of records in baseline, restraint stress, and shaker stress was 20 minutes, 60 minutes, and 20 minutes, respectively. Such a wide variety of data (different types of rats and different stressing surroundings) ensured the correct validation of the derived methods. Very slow signal component (trend), an outcome of free rat movements, was eliminated using a high-pass filter designed specifically for biomedical signals [26]. The cut-off frequency fH is tuned to the mean pulse interval of a subject: fH[Hz]=0.011/(mean(PI)[s]). Its value ranges from 0.055 Hz to 0.085 Hz in this experiment. The time series passed the stationarity, so the remaining signals were taken from 12 Wistar normotensive rats (total of 24 signals: 12 baseline, 6 restraint, and 6 shaker) and from 13 borderline hypertensive rats (total of 26 signals: 13 baseline, 7 restraint, and 6 shaker signal sets). All experimental procedures in this study conformed to the European Communities Council Directive of 24 November 1986 (86/609/ECC) and the Guidelines on Animal Experimentation of the School of Medicine, University of Belgrade.Short data series (number of samplesN≤750) were recorded from 41 medically examined male healthy volunteers positioned in a hospital bed, using Task Force® Monitor with sampling frequency equal to 1000 Hz (CNSystems Medizintechnik AG, Graz, http://www.cnsystems.com/products/task-force-monitor). Systolic blood pressure was obtained as local maxima of blood pressure waveforms, while R peaks were found as local maxima of electrocardiogram waveforms. The intervals between the successive R peaks form the R-R interval time series. Sensors were attached by medical professionals and patients were motionless, ensuring the reduced number of artifacts and stationarity of the signals. Short data lengths are typical for ambulatory monitoring, up to 10 minutes, and also ensure tranquility of the patient and stationarity of the recorded data. The experiment followed the ethical standards stated by the School of Medicine, University of Belgrade, and a signed permission was collected from each volunteer.Artificial data included samples with uniform distributions, generated using a reliable proprietary machine independent random number generator, and the samples with normal distribution generated using the polar method. For each one of the parameter sets, a hundred of the artificial data series were generated, each one comprising up toN=5000 data samples. Program code was written in Matlab 2014b.An illustration of the time series referenced within this study is shown in Figure1.Figure 1
Examples of real cardiovascular signals and of artificially generated signals. (a) SBP signal (gray) and PI signal (green) recorded from a rat at baseline condition, exposed to restraint stress and exposed to shaker stress. (b) SBP (gray) and R-R (green) signals recorded from a healthy volunteer. (c, d) Artificially generated signals with zero mean and standard deviation equal to 1; normal distribution (c) and uniform distribution (d).
(a)
(b)
(c)
(d)


### 2.2.XApEn in Brief

The procedure forXApEn estimation [22] implies xi∈X,i=1,…,N, as a “master” time series and yj∈Y,j=1,…,N, as its “follower,” where N is the number of signal samples. In biomedical data, SBP is considered as a master signal X, and PI is considered as its follower Y. The series are divided into the vectors of length m:(1)Xmi=xi,xi+τ,…,xi+m-1·τ,i=1,…,N-m-1·τYmj=yj,yj+τ,…,yj+m-1·τ,j=1,…,N-m-1·τ.Parameterτ can be introduced for vector decorrelation [21] but in most studies it is set to 1 [beat]. A vector from the master series is usually regarded as a “template.”Distance between two vectors is defined as a maximal absolute sample difference:(2)dXmi,Ymj=maxk=0,…,m-1⁡xi+k·τ-yj+k·τ,i,j=1,…,N-m-1·τ.The vectors are considered as “similar” if their distance is less than or equal to the predefined threshold (tolerance) r. The procedure of finding the similarities is called “template matching.” Standard scoring of time series must be performed; otherwise, vectors from different sources would be incomparable [22]. Standard scoring implies centralization and normalization of series x, (x-μ)/σ, where μ denotes mean and σ denotes standard deviation. For a correct estimation of μ and σ, the observed time series must be stationary at least in a wide sense.A conditional probability that a vectorY is within the distance r from a particular template Xm(i),i=1,…,N-(m-1)·τ, is estimated as follows:(3)p^imr=Pr⁡dXm,Ym≤r∣Xm=Xmi,Ym∈Y,r>0=1N-m-1·τ∑j=0N-m-1·τIdXmi,Ymj≤r,where “∧” denotes an estimate and I{·} denotes an indicator function that is equal to 1 if the condition it indicates is fulfilled (distance less than or equal to r). Denoting the number of “ones” with N1(i), the conditional probability (3) can be also expressed as(4)p^imr=N1iN-m-1·τ,N1i=∑j=0N-m-1·τIdXmi,Ymj≤r,i=1,…,N-m-1·τ.Probabilities (3) and (4) describe a binary event, as the vector Ym(j) can either match or mismatch its template Xm(i). A precise mathematical derivation of the number of binary events necessary for reliable probability estimation is given in [27]. The derivation in [27] is based assuming a normal approximation of binomial distribution, in our case observing the number of matches found during N-m+1 template matching trials. The resulting tables and graphs point out the probability that the estimated value p^im(r), for a given number of matching trials, lies within specific boundaries known as “confidence interval.” From the graphically presented results in [27], it can be seen that a 95% confidence interval of about (0.55·pim(r), 1.8·pim(r)), which is generally considered to be a reasonable uncertainty, is assessed if the number of trials is 10/pim(r). So, if the number of tests exceeds 10/pim(r), it is considered as a “weak criterion” that yields satisfactory good estimates of binary probabilities [27]. The number of matches that is equal to 100/pim(r) is considered as a “strong criterion” for which the boundaries are tighter, with probability estimation less uncertain, but usually it is more difficult to be achieved as it requires long time series. The corresponding confidence intervals, both for weak and for strong criteria, are derived in [27] as follows:(5)Pr⁡0.55·pimr≤p^imr≤1.8·pimr=0.95weak  criterion;Pr⁡0.8·pimr≤p^imr≤1.25·pimr=0.95strong  criterion.The next step inXApEn estimation is to average logarithms of the probability estimates (4) over all template values and to form a summand Φ^(m):(6)Φ^mr,N,τ=1N-m-1·τ∑i=0N-m-1·τln⁡p^imr.The procedure is repeated for the vectors of length m+1, yielding the estimate:(7)X^ApEnm,r,N,τ=Φ^mr,N,τ-Φ^m+1r,N,τ.Due to the fact that inXApEn two different time series are compared, a self-match (matching a template to itself), alleged source of bias inApEn, cannot occur. Unfortunately, comparison of vectors from different series also means that there is no guarantee that a template of length m in the master time series has at least one match in the follower time series; hence, a lack of matches p^im(r)=0 can occur. The consequence is a logarithm of zero in (6) that makesXApEn estimate impossible, pointed out in [5] as well. A usual correction ignores and skips p^im(r)=0 cases. Such a correction decreases the numerator in (6) with respect to its denominator and considerably underestimates the summand Φ^. We propose two summand Φ^(m) corrections, both different from the corrections in [5]. The first one is analogous toApEn and implements artificial self-matches p^im(r)=1/N-(m-1)·τ (see (8)). The second correction excludes p^im(r)=0 cases from the numerator of (6) (see (9)). If N0=∑i=0N-(m-1)·τI{p^im(r)=0} denotes a total number of zero matches, then(8)Φ^COR1mr,N,τ=∑i=0,p^imr≠0N-m-1·τln⁡p^imr-N0·ln⁡N-m-1·τN-m-1·τ(9)Φ^COR2mr,N,τ=∑i=0,p^imr≠0N-m-1·τln⁡p^imrN-m-1·τ-N0.A simple abstract experiment analyzes the proposed corrections. Suppose that all nonzero probabilities are equal to an arbitrary constant,p^im(r)=e-k,k>0. Then, summands Φ^ given by (6) with p^im(r)=0 skipped and given by (8) and by (9) are equal to(10)Φ^mr,N,τ=-k·1-N0N-m-1·τ6withp^imr=0skipped;-k·1-N0N-m-1·τ-N0·ln⁡N-m-1·τN-m-1·τ,correction8;-kcorrection9.If N tends to infinity, N0 converges to zero and all Φ^ estimates converge to the same true value, that is, to -k. But if N is limited, the underestimation of Φ^ in (6) (with p^im(r)=0 skipped) and (8) is proportional to N0. Equation (8) additionally suffers from a nonlinear bias due to the false self-matching, thus repeating the problem of genuine self-matching inApEn. The preferred correction is given by (9) and it was used in further evaluation. However, as it would be clarified later on, if theXApEn parameters are properly chosen, the need for this adjustment is reduced.


### 2.3. Parameter-InducedXApEn Inconsistency

The parameters necessary forXApEn estimation—time series length N, threshold level r, vector size m, and time lag τ—create a huge working space, but with unstable regions. RecentApEn contributions have reviewed this problem [4, 19], earlier also noticed in [6, 8, 15–18, 20]. But parameter study forXApEn does not exist, except for a brief comment in [5] (to the best of the authors’ knowledge). A few contributions that implementXApEn follow the standard guidelines from the originalApEn/XApEn papers [1–3, 22]: parameters m=1 and r=0.2 are implemented in [27]; parameters m=2 or 3, r=0.15, and N=1000 are implemented in [28, 29], while parameters m=1 and r=0.35 are used in [30]. It should be noted that the threshold value r inApEn studies is calculated with respect to the standard deviation σ, for example, r=0.2·σ. InXApEn studies, standard deviation is by definition set to σ=1, due to the obligatory standard scoring [22].An illustrative example of contradictory results induced by different working points is presented in Figure2. The contradictory results are noticed in borderline hypertensive rats exposed to shaker stress. Embedded plots show a “XApEn threshold profile” [17]:XApEn values estimated for a range of threshold r values, keeping all the other parameters fixed. To obtain a profile, threshold level was gradually increased with increments equal to 0.01 for r<0.2 and equal to 0.05 for r≥0.2. The black line corresponds toXApEn estimated from SBP and PI signals of rats in baseline conditions, while the red line corresponds to the threshold profile estimated from rats at stress. An intersection of these two lines divides the plot area into two regions: in the “gray” region,XApEn values are higher in baseline conditions than in stress, while in the “white” regionXApEn contradictorily becomes lower at baseline than at stress. The maximal value ofXApEn (the maximum of the threshold profile line) is denoted by a circle. An “empty” circle denotes the lower one of these two maximums. A complete threshold-length (r-N) plane is shown in the remaining part of Figure 2. A threshold profile is estimated for each particular time series length N. The thresholds positions rmax for whichXApEn profile reaches its maximum [17] are shown by circles. We have assumed a time lag τ=1.Figure 2
Illustrative example of inconsistent results: a threshold-length (r-N) plane forXApEn estimated from SBP and PI signals of rats at baseline (black lines and circles) and stress (red lines and circles). Gray regions: estimatedXApEn is higher at baseline than at stress. White regions:XApEn is higher at stress than at baseline. Circles: threshold rmax for whichXApEn reaches its maximal value; for a particular N, the empty circle denotes the lower of two maximalXApEn values. Embedded plot shows a complete threshold profile for N=4000. Results are obtained from borderline hypertensive rats exposed to shaker stress.The symbolic presentation in Figure2 distinguishes three different sources of contradictory results (flip-flop effect):(A)
For a fixed lengthN, theXApEn “flips” from being lower during stress for lower threshold r to being higher during stress for higher threshold r; for example, for N=1000, this change occurs for r≈0.15.(B)
For a fixed thresholdr, an undesirable change ofXApEn from being lower at stress to being higher at stress occurs at particular time series lengths; for example, if r=0.12, this change occurs between N=3000 and N=4000.(C)
Forrmax chosen for a threshold (rmax is a threshold for whichXApEn reaches its maximal value, [15–17]),XApEn changes from being higher at stress, if N=1000, shown by a full red circle, to being lower at stress for N>2000, shown by empty red circles.


### 2.4. Explicit Formulae for Error Estimation

The constitutiveXApEn elements are conditional probabilities p^im(r) ((3) and (4)), which are subject to an estimation error εim(r) each: (11)εimr=pimr-p^imr⟹p^imr=pimr·1-εimrpimr,i=1,…,N-m-1·τ.These errors are logarithmically accumulating to form a summand Φ^(m)(r,N,τ) estimation error εΦm(r,N,τ),(12)Φ^mr,N,τ=1N-m-1·τ∑i=0N-m-1·τln⁡pimr+ln⁡1-εimrpimr⟹εΦmr,N,τ=-1N-m-1·τ∑i=0N-m-1·τln⁡1-εimrpimr,andXApEn error εXm(r,N,τ),(13)εXmr,N,τ=1N-m·τ∑i=0N-m·τln⁡1-εim+1rpim+1r-1N-m-1·τ∑i=0N-m-1·τln⁡1-εimrpimr.To find a stable parameter region where estimation errors would be minimized, the true conditional probabilitiespim(r) must be evaluated for each template Xm(i): (14)pimr=∏k=0m-1∫xi+k·τ-rxi+k·τ+rfyy·dy,i=1,…,N-m-1·τ,where fy(y) denotes the probability density function (p.d.f.) of the signal Y (follower). A probability ∫xi+k·τ-rxi+k·τ+rfy(y)·dy corresponds to an event where a sample yj+k·τ is within the distance ±r from a template sample xi+k·τ.Unfortunately, biomedical signals are an outcome of numerous complex physiological inputs and feedback, so exact formulae for their p.d.f. have not been determined yet. Instead, we used artificially generated data with normal and uniform distribution:(15)fyUy=12·3,-3≤y≤30,elsewhere,uniform  p.d.f.fyNy=12·π·e-y2/2,normal  p.d.f.The artificial data were used both for the master and for the follower series, generating four experimental data sets: a master and its follower both normal, both uniform, normal master with uniform follower, and uniform master with normal follower.The corresponding true values of conditional probabilities, assumingτ=1, are(16)pimr=∏k=0m-1min⁡3,xi+k+r-max⁡-3,xi+k-r2·3m,xi+k-r<3,uniform  p.d.f.0,xi+k-r≥3uniform  p.d.f.∏k=0m-112·π∫xi+k-rxi+k+re-y2/2·dy=∏k=0m-1φxi+k+r-φxi+k-rnormal  p.d.f.Normal cumulative distribution function (1/2·π)∫-∞ze-y2/2·dy in the lower part of (16) is denoted by φ(z) instead of Φ(z), to avoid confusion with the summand Φ. With (16) known, all the relevant estimation errors (11)–(13) can be found.


### 2.5. Reference Thresholds for Cardiovascular Signals

The error analysis in artificial data has a privilege to evaluate a reference point (true value) using an exact formula. The problem with natural signals is that their data samples are statistically dependent. Besides, the signal nature is not completely revealed and there are no mathematical expressions to describe it. So, there exists no calculated reference value to which the estimated results could be compared. The first auxiliary point of reference consideringApEn was empirically obtained in [17], where the authors pointed out the existence ofApEn maximum and the corresponding threshold rmax for which this maximum is reached. rmax is found by performing a tedious task of estimatingApEn for each threshold value separately, thus obtaining a threshold profile. In [15], tedious estimation is substituted by a set of expressions that automatically evaluate a theoretical threshold rTH-A that is (almost) equal to rmax. In [16], the same authors proposed an alternative set of formulae for rTH-A. It was shown [15, 16] that the empirical threshold level rmax for which maximal value ofApEn is obtained when plotting the threshold profile depends on the time series length N, on standard deviations σ of time series x(i), and on standard deviations σD of lagged differential signal (x(i)-x(i-τ)), where τ is a time lag (delay) between the samples, τ=1. The observed dependency was shown to be linear with respect to the square root of σD and σ ratio, so the authors performed a computer search for the coefficients that yield the best fit in the least-squares sense to obtain (18). We have adopted the same method to find the m=1 case of a theoretical threshold rTH-A in (18), necessary for our evaluations, but not covered in [15, 16].A method for automatic evaluation of auxiliary reference point does not exist forXApEn. We observed [25] that the threshold rTH-A of master time series influences the threshold rmax ofXApEn and that their difference is linearly related to the square root of σDX+σDY, where σDX and σDY are standard deviations of the lagged differential series (x(i)-x(i-τ)) and (y(i)-y(i-τ)) [25]. Then, for each artificially created master X-follower Y combination (normal-normal, normal-uniform, uniform-normal, and uniform-uniform), we generated 100 normalized signal pairs, in total 400 signal pairs per parameter set. For each signal pair, we found the maximum of the threshold profile and the corresponding rmax, and we also estimated σDX and σDY. Then, the least-squares fitting is performed with respect to rmax, thus obtaining particular coefficients a(m), b(m), and c(m). The procedure is repeated for different N values, from 300 to 5000 with an increment equal to 100. The results were averaged for each particular m∈{1,2,3,4} giving the coefficient sets a(m), b(m), and c(m). Since several sets had almost the same averaged least-squares error, ten out of 25 baseline signals from rats were randomly selected for justification, and the coefficient set with minimal absolute error (averaged over all values of N and m) was selected for a threshold rTH-X, for whichXApEn approaches its maximal value:(17)rTH-Xm=rTH-Am+am+bm+cm·σDX+σDY/2N/10004,a,b,c=0,-0.015,0.03,m=1-0.02,0,0.023,m=20,-0.006,0.043,m=30,-0.11,0.13,m=4,where rTH-A is evaluated from the master time series according to(18)rTH-Am=em+fm·σDX/σXN/10004,e,f=-0.01,0.05,m=1-0.02,0.23,m=2-0.06,0.43,m=3-0.11,0.65,m=4.Here,σX is standard deviation of a time series X (inXApEn it is always equal to one) while σDX and σDY are standard deviations of the lagged differential series (x(i)-x(i-τ)) and (y(i)-y(i-τ)). The expression for m=1 in (18) is presented for the first time in this study, while the expressions for m=2,3,4 are taken from [15]. As it is shown in Table 3, the expressions from [15] yield a lower amount of relative error for cardiovascular data used in this study, so they are more appropriate than the similar expressions derived in [16]. The subset of thresholds (17) for m=2 and m=3 was presented in [25].The thresholdrTH-X is only halfway towards the reliable working point. As it would be shown, if threshold rTH-A is applied, the reliability criteria stated by N>10/pim(r) (or 100/pim(r)) are not satisfied. A different interpretation of these criteria states that at least 10 (or 100) events should be found, in order to have the estimated event probability fitted within the confidence bounds (5) with the certainty of 95%.Obviously, increasedm reduces the matching probability of the template and its follower, and thus a threshold increase had to be proportional to it. The coefficients in (19) and (20) are obtained using the same data sets as for (17). In this case, the conditional probabilities were estimated, and the fitting was performed among the coefficients that yield less than 5% of unreliable probabilities, over all values of m and over all data lengths:(19)rXWm=0.2·m+1+rTH-Xm,m≤4,N≤500,rXWm=0.17·m+1+rTH-Xm,m≤4,500<N≤2000,rXWm=0.15·m+1+rTH-Xm,m≤4,N>2000.In the case of strong reliability criteria, a relationship of threshold (20) with respect to vector length m is nonlinear:(20)rXSm=m+1+105N2·rTH-Xm+5-m10,m≤4.For the sake of clarity, the different threshold values are summarized in Table1.Table 1
Threshold abbreviations and definitions.
Abbreviation
Definition
r
⋅
σ
ApEn threshold, expressed in units of time series
r
XApEn threshold, dimensionless
r
m
a
x
Threshold obtained from threshold profile, for which (X)ApEn reaches its maximal value [17]
r
TH-A
Automatically calculatedApEn threshold using (18)
r
TH-X
Automatically calculatedXApEn threshold using (17)
r
X
W
Automatically calculatedXApEn threshold for which at least 95% of the estimated probabilities satisfy weak criterion.
r
X
S
Automatically calculatedXApEn threshold for which at least 95% of the estimated probabilities satisfy strong criterion.


### 2.6. Reliability Criteria in Statistically Dependent Data

The confidence bounds (5) hold for statistically independent binary events [27]. Template matching in artificial data satisfies this constraint: observing an m=1 template [xk], its successive matching tests |xk-yi|≤r and |xk-yi+1|≤r at the respective positions i and i+1 are independent, as an absolute difference of a constant xk and independent samples yi and yi+1 does not destroy independency. A similar line of reasoning holds if m>1. It follows that the estimated probabilities p^im(r) have a scaled binomial distribution [27], so the reliability criteria and confidence intervals (5) can be applied. If the data were statistically dependent, the only change is an increase of confidence bounds, without affecting the estimated values or the reliability criteria [27]. However, true confidence bounds in real biomedical data cannot be defined as true pim(r) is unknown, so statistical dependency is of no importance for real data. The reliability still can be tested in a form of N-m+1>(10or100)/p^im(r), showing whether the obtained value p^im(r) is estimated from a sufficient amount of data.


## 2.1. Experimental Data

The study used three types of data: artificial, animal, and human. The purpose of artificial data was to test the stability of estimated values, comparing them to the exact values calculated by the derived formulae, and also to develop a new set of threshold equations. The data recorded from animals were also used for method development and for its repeated testing, enabled by the record length. The data recorded from human volunteers were used to validate the derived methods in typical real applications where the records are short.The animal data include systolic blood pressure (SBP) signals and pulse interval (PI) signals recorded from two groups of rats: the first group comprised outbred male Wistar normotensive rats, while the second group comprised outbred male borderline hypertensive rats. Ten days before the recordings, radiotelemetric probes (TA11PA-C40, DSI, Transoma Medical) were implanted in the abdominal aorta. Blood pressure (BP) waveforms were digitally recorded and transmitted using Dataquest A.R.T. 4.0 equipment (DSI, St. Paul, MN, USA) with sampling frequencyfs = 1000 Hz. The systolic blood pressure (SBP) was determined at the receiver as waveform local maxima. The pulse interval was determined as an interval between the successive moments of maximal gradients ΔBP/ΔTmax, that is, the local maximal change (Δ) of blood pressure (BP) waveform per sampling interval ΔT. Pulse interval is an acceptable alternative to R-R interval of an electrocardiogram. During the experiments, the animals were exposed to two types of stress: shaker stress, with rats positioned on a platform shaking at 200 cycles/min, and restraint stress, with rats placed in a Plexiglas restrainer tube (ID 5.5 cm with pores) in the supine position, with a detailed description given in [24]. These data were already used for other research purposes [20, 24, 25]. The duration of records in baseline, restraint stress, and shaker stress was 20 minutes, 60 minutes, and 20 minutes, respectively. Such a wide variety of data (different types of rats and different stressing surroundings) ensured the correct validation of the derived methods. Very slow signal component (trend), an outcome of free rat movements, was eliminated using a high-pass filter designed specifically for biomedical signals [26]. The cut-off frequency fH is tuned to the mean pulse interval of a subject: fH[Hz]=0.011/(mean(PI)[s]). Its value ranges from 0.055 Hz to 0.085 Hz in this experiment. The time series passed the stationarity, so the remaining signals were taken from 12 Wistar normotensive rats (total of 24 signals: 12 baseline, 6 restraint, and 6 shaker) and from 13 borderline hypertensive rats (total of 26 signals: 13 baseline, 7 restraint, and 6 shaker signal sets). All experimental procedures in this study conformed to the European Communities Council Directive of 24 November 1986 (86/609/ECC) and the Guidelines on Animal Experimentation of the School of Medicine, University of Belgrade.Short data series (number of samplesN≤750) were recorded from 41 medically examined male healthy volunteers positioned in a hospital bed, using Task Force® Monitor with sampling frequency equal to 1000 Hz (CNSystems Medizintechnik AG, Graz, http://www.cnsystems.com/products/task-force-monitor). Systolic blood pressure was obtained as local maxima of blood pressure waveforms, while R peaks were found as local maxima of electrocardiogram waveforms. The intervals between the successive R peaks form the R-R interval time series. Sensors were attached by medical professionals and patients were motionless, ensuring the reduced number of artifacts and stationarity of the signals. Short data lengths are typical for ambulatory monitoring, up to 10 minutes, and also ensure tranquility of the patient and stationarity of the recorded data. The experiment followed the ethical standards stated by the School of Medicine, University of Belgrade, and a signed permission was collected from each volunteer.Artificial data included samples with uniform distributions, generated using a reliable proprietary machine independent random number generator, and the samples with normal distribution generated using the polar method. For each one of the parameter sets, a hundred of the artificial data series were generated, each one comprising up toN=5000 data samples. Program code was written in Matlab 2014b.An illustration of the time series referenced within this study is shown in Figure1.Figure 1
Examples of real cardiovascular signals and of artificially generated signals. (a) SBP signal (gray) and PI signal (green) recorded from a rat at baseline condition, exposed to restraint stress and exposed to shaker stress. (b) SBP (gray) and R-R (green) signals recorded from a healthy volunteer. (c, d) Artificially generated signals with zero mean and standard deviation equal to 1; normal distribution (c) and uniform distribution (d).
(a)
(b)
(c)
(d)


## 2.2.XApEn in Brief

The procedure forXApEn estimation [22] implies xi∈X,i=1,…,N, as a “master” time series and yj∈Y,j=1,…,N, as its “follower,” where N is the number of signal samples. In biomedical data, SBP is considered as a master signal X, and PI is considered as its follower Y. The series are divided into the vectors of length m:(1)Xmi=xi,xi+τ,…,xi+m-1·τ,i=1,…,N-m-1·τYmj=yj,yj+τ,…,yj+m-1·τ,j=1,…,N-m-1·τ.Parameterτ can be introduced for vector decorrelation [21] but in most studies it is set to 1 [beat]. A vector from the master series is usually regarded as a “template.”Distance between two vectors is defined as a maximal absolute sample difference:(2)dXmi,Ymj=maxk=0,…,m-1⁡xi+k·τ-yj+k·τ,i,j=1,…,N-m-1·τ.The vectors are considered as “similar” if their distance is less than or equal to the predefined threshold (tolerance) r. The procedure of finding the similarities is called “template matching.” Standard scoring of time series must be performed; otherwise, vectors from different sources would be incomparable [22]. Standard scoring implies centralization and normalization of series x, (x-μ)/σ, where μ denotes mean and σ denotes standard deviation. For a correct estimation of μ and σ, the observed time series must be stationary at least in a wide sense.A conditional probability that a vectorY is within the distance r from a particular template Xm(i),i=1,…,N-(m-1)·τ, is estimated as follows:(3)p^imr=Pr⁡dXm,Ym≤r∣Xm=Xmi,Ym∈Y,r>0=1N-m-1·τ∑j=0N-m-1·τIdXmi,Ymj≤r,where “∧” denotes an estimate and I{·} denotes an indicator function that is equal to 1 if the condition it indicates is fulfilled (distance less than or equal to r). Denoting the number of “ones” with N1(i), the conditional probability (3) can be also expressed as(4)p^imr=N1iN-m-1·τ,N1i=∑j=0N-m-1·τIdXmi,Ymj≤r,i=1,…,N-m-1·τ.Probabilities (3) and (4) describe a binary event, as the vector Ym(j) can either match or mismatch its template Xm(i). A precise mathematical derivation of the number of binary events necessary for reliable probability estimation is given in [27]. The derivation in [27] is based assuming a normal approximation of binomial distribution, in our case observing the number of matches found during N-m+1 template matching trials. The resulting tables and graphs point out the probability that the estimated value p^im(r), for a given number of matching trials, lies within specific boundaries known as “confidence interval.” From the graphically presented results in [27], it can be seen that a 95% confidence interval of about (0.55·pim(r), 1.8·pim(r)), which is generally considered to be a reasonable uncertainty, is assessed if the number of trials is 10/pim(r). So, if the number of tests exceeds 10/pim(r), it is considered as a “weak criterion” that yields satisfactory good estimates of binary probabilities [27]. The number of matches that is equal to 100/pim(r) is considered as a “strong criterion” for which the boundaries are tighter, with probability estimation less uncertain, but usually it is more difficult to be achieved as it requires long time series. The corresponding confidence intervals, both for weak and for strong criteria, are derived in [27] as follows:(5)Pr⁡0.55·pimr≤p^imr≤1.8·pimr=0.95weak  criterion;Pr⁡0.8·pimr≤p^imr≤1.25·pimr=0.95strong  criterion.The next step inXApEn estimation is to average logarithms of the probability estimates (4) over all template values and to form a summand Φ^(m):(6)Φ^mr,N,τ=1N-m-1·τ∑i=0N-m-1·τln⁡p^imr.The procedure is repeated for the vectors of length m+1, yielding the estimate:(7)X^ApEnm,r,N,τ=Φ^mr,N,τ-Φ^m+1r,N,τ.Due to the fact that inXApEn two different time series are compared, a self-match (matching a template to itself), alleged source of bias inApEn, cannot occur. Unfortunately, comparison of vectors from different series also means that there is no guarantee that a template of length m in the master time series has at least one match in the follower time series; hence, a lack of matches p^im(r)=0 can occur. The consequence is a logarithm of zero in (6) that makesXApEn estimate impossible, pointed out in [5] as well. A usual correction ignores and skips p^im(r)=0 cases. Such a correction decreases the numerator in (6) with respect to its denominator and considerably underestimates the summand Φ^. We propose two summand Φ^(m) corrections, both different from the corrections in [5]. The first one is analogous toApEn and implements artificial self-matches p^im(r)=1/N-(m-1)·τ (see (8)). The second correction excludes p^im(r)=0 cases from the numerator of (6) (see (9)). If N0=∑i=0N-(m-1)·τI{p^im(r)=0} denotes a total number of zero matches, then(8)Φ^COR1mr,N,τ=∑i=0,p^imr≠0N-m-1·τln⁡p^imr-N0·ln⁡N-m-1·τN-m-1·τ(9)Φ^COR2mr,N,τ=∑i=0,p^imr≠0N-m-1·τln⁡p^imrN-m-1·τ-N0.A simple abstract experiment analyzes the proposed corrections. Suppose that all nonzero probabilities are equal to an arbitrary constant,p^im(r)=e-k,k>0. Then, summands Φ^ given by (6) with p^im(r)=0 skipped and given by (8) and by (9) are equal to(10)Φ^mr,N,τ=-k·1-N0N-m-1·τ6withp^imr=0skipped;-k·1-N0N-m-1·τ-N0·ln⁡N-m-1·τN-m-1·τ,correction8;-kcorrection9.If N tends to infinity, N0 converges to zero and all Φ^ estimates converge to the same true value, that is, to -k. But if N is limited, the underestimation of Φ^ in (6) (with p^im(r)=0 skipped) and (8) is proportional to N0. Equation (8) additionally suffers from a nonlinear bias due to the false self-matching, thus repeating the problem of genuine self-matching inApEn. The preferred correction is given by (9) and it was used in further evaluation. However, as it would be clarified later on, if theXApEn parameters are properly chosen, the need for this adjustment is reduced.


## 2.3. Parameter-InducedXApEn Inconsistency

The parameters necessary forXApEn estimation—time series length N, threshold level r, vector size m, and time lag τ—create a huge working space, but with unstable regions. RecentApEn contributions have reviewed this problem [4, 19], earlier also noticed in [6, 8, 15–18, 20]. But parameter study forXApEn does not exist, except for a brief comment in [5] (to the best of the authors’ knowledge). A few contributions that implementXApEn follow the standard guidelines from the originalApEn/XApEn papers [1–3, 22]: parameters m=1 and r=0.2 are implemented in [27]; parameters m=2 or 3, r=0.15, and N=1000 are implemented in [28, 29], while parameters m=1 and r=0.35 are used in [30]. It should be noted that the threshold value r inApEn studies is calculated with respect to the standard deviation σ, for example, r=0.2·σ. InXApEn studies, standard deviation is by definition set to σ=1, due to the obligatory standard scoring [22].An illustrative example of contradictory results induced by different working points is presented in Figure2. The contradictory results are noticed in borderline hypertensive rats exposed to shaker stress. Embedded plots show a “XApEn threshold profile” [17]:XApEn values estimated for a range of threshold r values, keeping all the other parameters fixed. To obtain a profile, threshold level was gradually increased with increments equal to 0.01 for r<0.2 and equal to 0.05 for r≥0.2. The black line corresponds toXApEn estimated from SBP and PI signals of rats in baseline conditions, while the red line corresponds to the threshold profile estimated from rats at stress. An intersection of these two lines divides the plot area into two regions: in the “gray” region,XApEn values are higher in baseline conditions than in stress, while in the “white” regionXApEn contradictorily becomes lower at baseline than at stress. The maximal value ofXApEn (the maximum of the threshold profile line) is denoted by a circle. An “empty” circle denotes the lower one of these two maximums. A complete threshold-length (r-N) plane is shown in the remaining part of Figure 2. A threshold profile is estimated for each particular time series length N. The thresholds positions rmax for whichXApEn profile reaches its maximum [17] are shown by circles. We have assumed a time lag τ=1.Figure 2
Illustrative example of inconsistent results: a threshold-length (r-N) plane forXApEn estimated from SBP and PI signals of rats at baseline (black lines and circles) and stress (red lines and circles). Gray regions: estimatedXApEn is higher at baseline than at stress. White regions:XApEn is higher at stress than at baseline. Circles: threshold rmax for whichXApEn reaches its maximal value; for a particular N, the empty circle denotes the lower of two maximalXApEn values. Embedded plot shows a complete threshold profile for N=4000. Results are obtained from borderline hypertensive rats exposed to shaker stress.The symbolic presentation in Figure2 distinguishes three different sources of contradictory results (flip-flop effect):(A)
For a fixed lengthN, theXApEn “flips” from being lower during stress for lower threshold r to being higher during stress for higher threshold r; for example, for N=1000, this change occurs for r≈0.15.(B)
For a fixed thresholdr, an undesirable change ofXApEn from being lower at stress to being higher at stress occurs at particular time series lengths; for example, if r=0.12, this change occurs between N=3000 and N=4000.(C)
Forrmax chosen for a threshold (rmax is a threshold for whichXApEn reaches its maximal value, [15–17]),XApEn changes from being higher at stress, if N=1000, shown by a full red circle, to being lower at stress for N>2000, shown by empty red circles.


## 2.4. Explicit Formulae for Error Estimation

The constitutiveXApEn elements are conditional probabilities p^im(r) ((3) and (4)), which are subject to an estimation error εim(r) each: (11)εimr=pimr-p^imr⟹p^imr=pimr·1-εimrpimr,i=1,…,N-m-1·τ.These errors are logarithmically accumulating to form a summand Φ^(m)(r,N,τ) estimation error εΦm(r,N,τ),(12)Φ^mr,N,τ=1N-m-1·τ∑i=0N-m-1·τln⁡pimr+ln⁡1-εimrpimr⟹εΦmr,N,τ=-1N-m-1·τ∑i=0N-m-1·τln⁡1-εimrpimr,andXApEn error εXm(r,N,τ),(13)εXmr,N,τ=1N-m·τ∑i=0N-m·τln⁡1-εim+1rpim+1r-1N-m-1·τ∑i=0N-m-1·τln⁡1-εimrpimr.To find a stable parameter region where estimation errors would be minimized, the true conditional probabilitiespim(r) must be evaluated for each template Xm(i): (14)pimr=∏k=0m-1∫xi+k·τ-rxi+k·τ+rfyy·dy,i=1,…,N-m-1·τ,where fy(y) denotes the probability density function (p.d.f.) of the signal Y (follower). A probability ∫xi+k·τ-rxi+k·τ+rfy(y)·dy corresponds to an event where a sample yj+k·τ is within the distance ±r from a template sample xi+k·τ.Unfortunately, biomedical signals are an outcome of numerous complex physiological inputs and feedback, so exact formulae for their p.d.f. have not been determined yet. Instead, we used artificially generated data with normal and uniform distribution:(15)fyUy=12·3,-3≤y≤30,elsewhere,uniform  p.d.f.fyNy=12·π·e-y2/2,normal  p.d.f.The artificial data were used both for the master and for the follower series, generating four experimental data sets: a master and its follower both normal, both uniform, normal master with uniform follower, and uniform master with normal follower.The corresponding true values of conditional probabilities, assumingτ=1, are(16)pimr=∏k=0m-1min⁡3,xi+k+r-max⁡-3,xi+k-r2·3m,xi+k-r<3,uniform  p.d.f.0,xi+k-r≥3uniform  p.d.f.∏k=0m-112·π∫xi+k-rxi+k+re-y2/2·dy=∏k=0m-1φxi+k+r-φxi+k-rnormal  p.d.f.Normal cumulative distribution function (1/2·π)∫-∞ze-y2/2·dy in the lower part of (16) is denoted by φ(z) instead of Φ(z), to avoid confusion with the summand Φ. With (16) known, all the relevant estimation errors (11)–(13) can be found.


## 2.5. Reference Thresholds for Cardiovascular Signals

The error analysis in artificial data has a privilege to evaluate a reference point (true value) using an exact formula. The problem with natural signals is that their data samples are statistically dependent. Besides, the signal nature is not completely revealed and there are no mathematical expressions to describe it. So, there exists no calculated reference value to which the estimated results could be compared. The first auxiliary point of reference consideringApEn was empirically obtained in [17], where the authors pointed out the existence ofApEn maximum and the corresponding threshold rmax for which this maximum is reached. rmax is found by performing a tedious task of estimatingApEn for each threshold value separately, thus obtaining a threshold profile. In [15], tedious estimation is substituted by a set of expressions that automatically evaluate a theoretical threshold rTH-A that is (almost) equal to rmax. In [16], the same authors proposed an alternative set of formulae for rTH-A. It was shown [15, 16] that the empirical threshold level rmax for which maximal value ofApEn is obtained when plotting the threshold profile depends on the time series length N, on standard deviations σ of time series x(i), and on standard deviations σD of lagged differential signal (x(i)-x(i-τ)), where τ is a time lag (delay) between the samples, τ=1. The observed dependency was shown to be linear with respect to the square root of σD and σ ratio, so the authors performed a computer search for the coefficients that yield the best fit in the least-squares sense to obtain (18). We have adopted the same method to find the m=1 case of a theoretical threshold rTH-A in (18), necessary for our evaluations, but not covered in [15, 16].A method for automatic evaluation of auxiliary reference point does not exist forXApEn. We observed [25] that the threshold rTH-A of master time series influences the threshold rmax ofXApEn and that their difference is linearly related to the square root of σDX+σDY, where σDX and σDY are standard deviations of the lagged differential series (x(i)-x(i-τ)) and (y(i)-y(i-τ)) [25]. Then, for each artificially created master X-follower Y combination (normal-normal, normal-uniform, uniform-normal, and uniform-uniform), we generated 100 normalized signal pairs, in total 400 signal pairs per parameter set. For each signal pair, we found the maximum of the threshold profile and the corresponding rmax, and we also estimated σDX and σDY. Then, the least-squares fitting is performed with respect to rmax, thus obtaining particular coefficients a(m), b(m), and c(m). The procedure is repeated for different N values, from 300 to 5000 with an increment equal to 100. The results were averaged for each particular m∈{1,2,3,4} giving the coefficient sets a(m), b(m), and c(m). Since several sets had almost the same averaged least-squares error, ten out of 25 baseline signals from rats were randomly selected for justification, and the coefficient set with minimal absolute error (averaged over all values of N and m) was selected for a threshold rTH-X, for whichXApEn approaches its maximal value:(17)rTH-Xm=rTH-Am+am+bm+cm·σDX+σDY/2N/10004,a,b,c=0,-0.015,0.03,m=1-0.02,0,0.023,m=20,-0.006,0.043,m=30,-0.11,0.13,m=4,where rTH-A is evaluated from the master time series according to(18)rTH-Am=em+fm·σDX/σXN/10004,e,f=-0.01,0.05,m=1-0.02,0.23,m=2-0.06,0.43,m=3-0.11,0.65,m=4.Here,σX is standard deviation of a time series X (inXApEn it is always equal to one) while σDX and σDY are standard deviations of the lagged differential series (x(i)-x(i-τ)) and (y(i)-y(i-τ)). The expression for m=1 in (18) is presented for the first time in this study, while the expressions for m=2,3,4 are taken from [15]. As it is shown in Table 3, the expressions from [15] yield a lower amount of relative error for cardiovascular data used in this study, so they are more appropriate than the similar expressions derived in [16]. The subset of thresholds (17) for m=2 and m=3 was presented in [25].The thresholdrTH-X is only halfway towards the reliable working point. As it would be shown, if threshold rTH-A is applied, the reliability criteria stated by N>10/pim(r) (or 100/pim(r)) are not satisfied. A different interpretation of these criteria states that at least 10 (or 100) events should be found, in order to have the estimated event probability fitted within the confidence bounds (5) with the certainty of 95%.Obviously, increasedm reduces the matching probability of the template and its follower, and thus a threshold increase had to be proportional to it. The coefficients in (19) and (20) are obtained using the same data sets as for (17). In this case, the conditional probabilities were estimated, and the fitting was performed among the coefficients that yield less than 5% of unreliable probabilities, over all values of m and over all data lengths:(19)rXWm=0.2·m+1+rTH-Xm,m≤4,N≤500,rXWm=0.17·m+1+rTH-Xm,m≤4,500<N≤2000,rXWm=0.15·m+1+rTH-Xm,m≤4,N>2000.In the case of strong reliability criteria, a relationship of threshold (20) with respect to vector length m is nonlinear:(20)rXSm=m+1+105N2·rTH-Xm+5-m10,m≤4.For the sake of clarity, the different threshold values are summarized in Table1.Table 1
Threshold abbreviations and definitions.
Abbreviation
Definition
r
⋅
σ
ApEn threshold, expressed in units of time series
r
XApEn threshold, dimensionless
r
m
a
x
Threshold obtained from threshold profile, for which (X)ApEn reaches its maximal value [17]
r
TH-A
Automatically calculatedApEn threshold using (18)
r
TH-X
Automatically calculatedXApEn threshold using (17)
r
X
W
Automatically calculatedXApEn threshold for which at least 95% of the estimated probabilities satisfy weak criterion.
r
X
S
Automatically calculatedXApEn threshold for which at least 95% of the estimated probabilities satisfy strong criterion.


## 2.6. Reliability Criteria in Statistically Dependent Data

The confidence bounds (5) hold for statistically independent binary events [27]. Template matching in artificial data satisfies this constraint: observing an m=1 template [xk], its successive matching tests |xk-yi|≤r and |xk-yi+1|≤r at the respective positions i and i+1 are independent, as an absolute difference of a constant xk and independent samples yi and yi+1 does not destroy independency. A similar line of reasoning holds if m>1. It follows that the estimated probabilities p^im(r) have a scaled binomial distribution [27], so the reliability criteria and confidence intervals (5) can be applied. If the data were statistically dependent, the only change is an increase of confidence bounds, without affecting the estimated values or the reliability criteria [27]. However, true confidence bounds in real biomedical data cannot be defined as true pim(r) is unknown, so statistical dependency is of no importance for real data. The reliability still can be tested in a form of N-m+1>(10or100)/p^im(r), showing whether the obtained value p^im(r) is estimated from a sufficient amount of data.


## 3. Results and Discussion

### 3.1. Estimation Errors in Artificial Data

ForXApEn error estimation, master series X and a hundred of its followers Y were artificially generated for a range of (r,N,m,(τ=1)) parameters, and the same values exactly were calculated using the derived formulae. The results for all four data sets, due to the size, are presented within the Supplementary Material (available online at https://doi.org/10.1155/2017/8365685).Relative conditional probability errorsεim(r)/pim(r)[%] are plotted in Figure 3, for series X and Y both normally distributed and for τ=1. Satisfactory results are obtained for an extremely atypical threshold level r=1.50. For a typical threshold r=0.15, error is below 10% only for large N and m=1. The results are in accordance with reliability criteria N-m+1>(10or100)/pim(r) stating that N or pim(r) should be large. Indeed, pim(r) is defined as a product of m probabilities (14), so it is at its largest if m=1, while increased threshold r augments the number of template matches N1(i) (4), again enlarging pim(r).Figure 3
The relative conditional probability error[%] in artificial data for a typical threshold value r=0.15 (a), for a higher threshold r=0.60 (b), and for an unusually high threshold r=1.50 (c); the gray dashed horizontal line marks 10% error; results are presented in a log-log plane as mean±σ (standard deviation). Artificial master and follower series are both normally distributed.
(a)
(b)
(c)The reliability ofp^im(r) estimate is checked comparing the data length N-m+1 with the required minimum 10/pim(r). The percentage of test failures is presented in Figure 4, in an r-N plane, where 0% failures are presented by green color that gradually changes to yellow (100% failures). The white line marks classical r and N values, revealing that, for m=3 and 4, the estimates that implement classical parameters are unreliable. Nonlinear interrelation of r, N, and m in Figure 4 also points out that the choice of parameters that would yield a stable working point is not a straightforward task.Figure 4
Percentage of estimatesp^im(r) in artificial data that fail to satisfy the criterion N-m+1≥10/p^im(r)[%], presented in threshold-length (r-N) plane, for m=1,2,3,4; failures are presented by color changes from green (0% failures) to yellow (100% failures); (a) X and Y normally distributed; (b) X and Y uniformly distributed; white horizontal and vertical lines cut each plane at typical parameter values r=0.2 and N=1000.
(a)
(b)The estimation error of summandΦ^ (12) exposes two problems. The first one is a lack of matches (p^im(r)=0) that induce logarithms of zero in (6) which is, as shown in Figure 5, a subset of all unreliable probabilities from Figure 4. The second problem is an error accumulation (12).Figure 5
Number ofp^im(r)=0 estimates [%] in artificial data presented in threshold-length (r-N) plane for m=3 and 4; color changes from green (0% zero matches) to yellow (100% zero matches); (a) X and Y normally distributed; (b) X and Y uniformly distributed.
(a)
(b)Figure6 shows a line graph of summand Φ and its changes with respect to series length N if it is evaluated by a formula (exact values), if it is estimated in a traditional way, and if a correction (9) is applied. This correction does attenuate estimation errors, even for high m values, but it cannot be used alone, as it removes a subset of the problem only. The convergence of summand Φ^ towards its true value depends on data type: uniform data achieve the correct limit for shorter lengths N. The convergence is also m-dependent, requiring different data lengths to achieve the same limit. Since X^ApEn is defined as a difference of Φ^ values (7), the fact that each summand Φ^ adds its own and independent estimation error is a further source of X^ApEn instability.Figure 6
Summand Φ forr=0.3: estimated omitting the zero probability, estimated and corrected, and calculated exactly using a formula; (a) X and Y normally distributed; (b) X and Y uniformly distributed.
(a)
(b)Figure7 shows two-dimensionalXApEn profiles. Values calculated using the derived formulae ((a) and gray lines in (b) and (c)) verify the known fact thatXApEn linearly decreases with the logarithm of threshold value [5]. The major conclusion is that the trueXApEn values are influenced neither by time series length N nor by template length m. Frequent notions thatXApEn monotonously increases with N until a plateau is reached or that its value decreases with m are a consequence of probability underestimation that decreases with series length and eventually disappears. The panels in the second row of (a) show corrected X^ApEn estimates that are in good accordance with true values for m=1, but showing increased discrepancies for m=2 and 3. If no correction is applied (the third row of (a)),    X^ApEn substantially diverges from the correct values. (b) and (c) show X^ApEn as a line graph. The discrepancy between the true and estimated values points out that generally accepted guidelines for parameter choice do not offer a safe X^ApEn estimation.Figure 7
XApEn of artificial data with normal distribution in threshold-length (r-N) plane, for m=1,2, and 3. The three rows of (a) show, respectively, calculated (true)XApEn,XApEn estimated with correction, andXApEn estimated omitting log⁡(0). Dark green color corresponds toXApEn estimated within the unstableN-r-m parameter region. Light green indicates lowXApEn, while yellow color indicates highXApEn. Horizontal cuts at N=2000 are presented as line graph in (b); vertical cuts at r=0.3 are presented as line graph in (c).
(a)
(b)
(c)


### 3.2. Stable Working Point in Cardiovascular Environment

To verify the accuracy of the derived equation (17), we followed the recommendations from [15, 16], where their equations are validated with respect toApEn relative error. For each one of SBP-PI (or R-R) and PI (or R-R)-SBP signal pairs, of animals (or volunteers), a threshold rTH-X was evaluated according to (17) and used as a parameter forXApEn(rTH-X) estimation. The true maximalXApEn(rmax) value is provided by a time-consuming procedure of finding a threshold profile. The corresponding relative error εXA_max is presented in Table 2. The results estimated from signals recorded in human volunteers are presented in the italic part of Table 2. The relative error is below 2%, except for very short data lengths (N=300) and m=1 when it reaches 6%; unfortunately, this corresponds to the ambulatory records that cannot be too long.Table 2
Relative error ofXApEn(rTH-X) with respect to XApEn(rmax) in real signals [%].
N
εXA_max=XApEn(rTH-X)-XApEn(rmax)/XApEn(rmax)⋅100
X
=
P
I, Y=SBP
X
=
S
B
P, Y=PI
m = 1
m = 2
m = 3
m = 4
m = 1
m = 2
m = 3
m = 4
300
5.16 ± 0.40
1.25 ± 0.28
1.61 ± 0.21
5.51 ± 0.75
6.40 ± 0.61
2.42 ± 0.49
2.85 ± 0.65
5.60 ± 0.88
500
2.57 ± 0.25
1.13 ± 0.21
2.61 ± 0.51
2.91 ± 0.57
3.61 ± 0.39
1.74 ± 0.29
1.68 ± 0.34
2.83 ± 0.60
750
1.57 ± 0.18
0.90 ± 0.21
1.93 ± 0.49
2.06 ± 0.44
2.09 ± 0.26
1.09 ± 0.25
1.56 ± 0.28
1.62 ± 0.43
2000
0.64 ± 0.17
1.07 ± 0.26
1.50 ± 0.45
1.83 ± 0.54
0.49 ± 0.16
1.43 ± 0.29
1.96 ± 1.40
2.62 ± 0.47
2500
0.46 ± 0.07
0.88 ± 0.14
1.29 ± 0.29
1.61 ± 0.42
0.32 ± 0.05
1.01 ± 0.22
1.45 ± 0.29
2.99 ± 0.47
3000
0.80 ± 0.16
0.91 ± 0.18
0.99 ± 0.25
1.22 ± 0.26
0.33 ± 0.08
0.99 ± 0.18
0.94 ± 0.16
2.03 ± 0.38
3500
0.67 ± 0.12
1.04 ± 0.19
1.89 ± 0.67
2.68 ± 1.28
0.52 ± 0.14
1.41 ± 0.33
1.58 ± 0.47
2.64 ± 0.58
4000
1.03 ± 0.20
1.65 ± 0.28
1.58 ± 0.26
1.15 ± 0.23
1.08 ± 0.25
1.22 ± 0.23
1.09 ± 0.24
1.79 ± 0.45
4500
0.68 ± 0.18
1.41 ± 0.39
1.72 ± 0.59
2.80 ± 0.73
0.94 ± 0.18
1.20 ± 0.30
1.90 ± 0.92
2.28 ± 1.33
5000
1.59 ± 0.36
1.57 ± 0.44
1.73 ± 0.61
2.26 ± 0.75
0.42 ± 0.08
2.20 ± 0.39
1.74 ± 0.56
1.64 ± 0.40
Results are presented as mean ± SEM. Results of 41 human volunteers are shown in the first three italic rows. Results of 25 rats in different conditions (total of 50 signals) are presented in roman rows.Table 3
Relative error ofApEn(rTH-A) with respect to ApEn(rmax) [%] in real signals.
N
εA_max=ApEn(rTH-A)-ApEn(rmax)/ApEn(rmax)⋅100
X
=
P
I
X
=
S
B
P
r
TH-A [15]
r
TH-A [16]
r
TH-A [15]
r
TH-A [16]
2000
3.46 ± 0.53
4.07 ± 0.54
1.95 ± 0.33
2.33 ± 0.35
3000
3.23 ± 0.54
3.82 ± 0.56
2.30 ± 0.35
2.91 ± 0.48
4000
3.72 ± 0.65
4.40 ± 0.67
2.86 ± 0.48
3.33 ± 0.51
5000
3.75 ± 0.67
6.45 ± 1.26
2.68 ± 0.43
5.63 ± 1.29
Error is evaluated according to the guidelines from [15] and from [16]. Results are presented as mean ± SEM, averaged over 100 time series: 50 data series of Wistar and borderline hypertensive rats in baseline and in stress, repeated for m=2 and for m=3.For the sake of comparison, we estimatedApEn error from our data as well. We applied rTH-A calculated according to [15] (see (18)) and according to [16]. The rTH-A relative error is shown in Table 3, providing the following conclusions: (1) error level is lower if the equations from [15] are used and (2) comparison of results in Tables 2 and 3 verifies that our expressions (17) are equal in performance to expressions (18) derived by Lu et al. in [15, 16].The valuerTH-X cannot be applied, since the percentage of the probabilities that fail the weak and strong criteria is almost 100% (Figure 8). But this threshold reflects the nature of data and it can be used as a reference for thresholds rXW and rXS. Their implementation reduces the failures, as shown in Figure 8.Figure 8
Unreliable estimated probabilitiesp^im(r)[%] in real signals; (a) r∈{rXW,rTH-X} and weak criterion is applied; (b) r∈{rXS,rTH-X} and strong criterion is applied. The gray dashed horizontal line marks 10% of the unreliable probability estimates. The left part of each panel shows the results of Wistar and borderline hypertensive rats at baseline and stress; the right part of each panel shows the results of volunteers.
(a)
(b)The actual threshold values calculated from real signals are shown in Table4. For the sake of comparison, threshold values for normally distributed artificial data are also plotted (Figure 9). These thresholds show that the traditional guidelines (0.15, 0.2, and 0.35) considerably underestimate the actual requirements for the reliableXApEn estimation. Even if m=1, the threshold required for strong reliability may exceed r=3. Although it was speculated [1] that detailed system information would be lost if r values were larger than 0.1·σ or 0.2·σ (forApEn), meaning that if the threshold exceeds these values a template would always match its follower, in investigated cardiovascular signals, this assumption does not hold as the evaluated high threshold levels are considerably below the amplitude range of the observed data.Table 4
Threshold values obtained from real signals.
(a)
ThresholdrTH-X for which XApEn reaches its maximal value (17)
N
X
=
P
I, Y=SBP
X
=
S
B
P, Y=PI
m = 1
m = 2
m = 3
m = 4
m = 1
m = 2
m = 3
m = 4
300
0.05 ± 0.01
0.24 ± 0.01
0.44 ± 0.01
0.62 ± 0.01
0.05 ± 0.00
0.22 ± 0.01
0.40 ± 0.01
0.55 ± 0.02
500
0.05 ± 0.01
0.21 ± 0.00
0.39 ± 0.01
0.53 ± 0.00
0.04 ± 0.00
0.19 ± 0.00
0.35 ± 0.01
0.47 ± 0.01
750
0.04 ± 0.01
0.18 ± 0.04
0.35 ± 0.01
0.47 ± 0.01
0.04 ± 0.01
0.16 ± 0.00
0.31 ± 0.01
0.41 ± 0.01
2000
0.05 ± 0.01
0.19 ± 0.00
0.37 ± 0.01
0.51 ± 0.01
0.04 ± 0.00
0.14 ± 0.00
0.27 ± 0.00
0.60 ± 0.01
2500
0.05 ± 0.01
0.18 ± 0.00
0.35 ± 0.01
0.48 ± 0.01
0.04 ± 0.00
0.13 ± 0.00
0.25 ± 0.01
0.34 ± 0.01
3000
0.04 ± 0.01
0.18 ± 0.00
0.33 ± 0.01
0.46 ± 0.01
0.03 ± 0.00
0.13 ± 0.00
0.24 ± 0.00
0.33 ± 0.05
3500
0.04 ± 0.01
0.17 ± 0.00
0.32 ± 0.01
0.45 ± 0.01
0.03 ± 0.00
0.12 ± 0.00
0.23 ± 0.00
0.31 ± 0.05
4000
0.04 ± 0.01
0.17 ± 0.00
0.31 ± 0.01
0.44 ± 0.01
0.03 ± 0.00
0.17 ± 0.00
0.22 ± 0.00
0.30 ± 0.05
4500
0.04 ± 0.01
0.11 ± 0.00
0.30 ± 0.01
0.43 ± 0.01
0.03 ± 0.00
0.11 ± 0.00
0.21 ± 0.00
0.29 ± 0.01
5000
0.04 ± 0.01
0.16 ± 0.00
0.29 ± 0.01
0.42 ± 0.01
0.03 ± 0.00
0.11 ± 0.00
0.20 ± 0.00
0.29 ± 0.01
(b)
ThresholdrXW for which at least 95% of the probabilities satisfy N-m+1>10/p^im(rXW) (19)
N
X
=
P
I, Y=SBP
X
=
S
B
P, Y=PI
m = 1
m = 2
m = 3
m = 4
m = 1
m = 2
m = 3
m = 4
300
0.46 ± 0.00
0.84 ± 0.01
1.24 ± 0.01
1.62 ± 0.01
0.45 ± 0.00
0.82 ± 0.01
1.20 ± 0.01
1.55 ± 0.02
500
0.45 ± 0.00
0.81 ± 0.00
1.19 ± 0.01
1.53 ± 0.01
0.44 ± 0.00
0.78 ± 0.00
1.15 ± 0.01
1.47 ± 0.01
750
0.38 ± 0.00
0.69 ± 0.00
1.03 ± 0.01
1.32 ± 0.01
0.38 ± 0.00
0.67 ± 0.00
0.99 ± 0.01
1.26 ± 0.01
2000
0.39 ± 0.00
0.70 ± 0.00
1.05 ± 0.01
1.36 ± 0.01
0.37 ± 0.00
0.65 ± 0.00
0.95 ± 0.01
1.21 ± 0.01
2500
0.35 ± 0.00
0.63 ± 0.00
0.95 ± 0.01
1.23 ± 0.01
0.34 ± 0.00
0.58 ± 0.00
0.85 ± 0.01
1.09 ± 0.01
3000
0.34 ± 0.00
0.63 ± 0.00
0.93 ± 0.01
1.22 ± 0.01
0.33 ± 0.00
0.59 ± 0.00
0.84 ± 0.00
1.08 ± 0.01
3500
0.34 ± 0.00
0.62 ± 0.00
0.92 ± 0.01
1.20 ± 0.01
0.33 ± 0.00
0.57 ± 0.00
0.83 ± 0.00
1.06 ± 0.01
4000
0.34 ± 0.00
0.62 ± 0.00
0.91 ± 0.01
1.19 ± 0.01
0.33 ± 0.00
0.57 ± 0.00
0.82 ± 0.00
1.05 ± 0.01
4500
0.34 ± 0.00
0.61 ± 0.00
0.90 ± 0.01
1.18 ± 0.01
0.33 ± 0.00
0.56 ± 0.00
0.81 ± 0.00
1.04 ± 0.01
5000
0.34 ± 0.00
0.61 ± 0.00
0.89 ± 0.01
1.17 ± 0.01
0.33 ± 0.00
0.56 ± 0.00
0.80 ± 0.00
1.04 ± 0.01
(c)
ThresholdrXS for which at least 95% of the probabilities satisfy N-m+1>100/p^im(rXS) (20)
N
X
=
P
I, Y=SBP
X
=
S
B
P, Y=PI
m = 1
m = 2
m = 3
m = 4
m = 1
m = 2
m = 3
m = 4
300
1.42 ± 0.00
2.21 ± 0.02
3.28 ± 0.04
4.37 ± 0.09
1.40 ± 0.01
2.12 ± 0.02
3.07 ± 0.05
3.99 ± 0.10
500
1.08 ± 0.00
1.73 ± 0.14
2.60 ± 0.04
3.41 ± 0.01
1.06 ± 0.00
1.65 ± 0.01
2.41 ± 0.00
3.06 ± 0.07
750
0.96 ± 0.00
1.54 ± 0.01
2.28 ± 0.03
2.93 ± 0.05
0.96 ± 0.00
1.47 ± 0.01
2.13 ± 0.03
2.64 ± 0.04
2000
0.91 ± 0.00
1.50 ± 0.01
2.29 ± 0.03
3.06 ± 0.05
0.88 ± 0.00
1.33 ± 0.01
1.89 ± 0.02
2.29 ± 0.00
2500
0.90 ± 0.00
1.46 ± 0.01
2.21 ± 0.03
2.93 ± 0.05
0.88 ± 0.00
1.31 ± 0.01
1.82 ± 0.02
2.20 ± 0.03
3000
0.89 ± 0.00
1.44 ± 0.00
2.14 ± 0.03
2.83 ± 0.04
0.87 ± 0.00
1.29 ± 0.01
1.77 ± 0.02
2.14 ± 0.03
3500
0.89 ± 0.00
1.41 ± 0.01
2.08 ± 0.03
2.76 ± 0.04
0.87 ± 0.00
1.27 ± 0.01
1.72 ± 0.02
2.07 ± 0.02
4000
0.88 ± 0.00
1.40 ± 0.01
2.05 ± 0.02
2.72 ± 0.04
0.86 ± 0.00
1.25 ± 0.01
1.67 ± 0.02
2.00 ± 0.03
4500
0.88 ± 0.00
1.39 ± 0.01
2.00 ± 0.02
2.66 ± 0.04
0.86 ± 0.00
1.24 ± 0.01
1.63 ± 0.02
1.95 ± 0.02
5000
0.88 ± 0.00
1.37 ± 0.01
1.97 ± 0.02
2.62 ± 0.04
0.86 ± 0.00
1.23 ± 0.01
1.60 ± 0.00
1.93 ± 0.02
Results are presented as mean ± SEM, averaged over 50 time series of both Wistar and borderline hypertensive rats, both in baseline and in stress, and averaged over 41 time series of volunteers; results from human volunteers are in italics.Figure 9
Calculated threshold values for artificial data, with signalsX and Y both normally distributed. (a) Threshold for whichXApEn gets maximal value (see (17)). (b) Threshold for which at least 95% of the probabilities satisfy weak criterion (see (19)). (c) Threshold for which at least 95% of the probabilities satisfy strong criterion (see (20)). The gray horizontal dashed line corresponds to r=0.3.
(a)
(b)
(c)The developed sets of threshold expressions are primarily intended for cardiovascular signal analysis. The artificial signals used as well are uniform on segment (with bounded amplitudes) and Gaussian (with increased incidence of small amplitudes, but unbounded). Both artificial distributions are symmetric. The amplitudes of real signals (PI, SBP, and R-R) are bounded by physiological constraints. Signal asymmetry is expressed by skewness factor, shown in Table5. Cardiovascular data can be skewed to the left and to the right, so Table 5 presents the average of absolute skewness values. SBP and PI signals exhibit a level of asymmetry comparable to Rayleigh distribution (skewness ≈ 0.631), but below the skewness of exponential distribution (skewness ≈ 2). Threshold evaluation in signals with highly skewed distribution would be a subject for follow-up studies. In particular, it would be precious to explore an important class of interspike interval time series, with distributions that are close to exponential. A memoryless property of exponential distribution might affect the values of differential time series standard deviation σD, necessary for (17) and (18).Table 5
Level of asymmetry of real signals (skewness).
Rats
Volunteers
SBP
PI
SBP
R-R
Absolute skewness
0.64 ± 0.86
0.48 ± 1.02
0.42 ± 0.54
0.35 ± 0.33
Portion of signals with negative skewness[%]
36%
24%
46%
60%
Absolute skewness is expressed asmean±σ; each signal set comprises a couple of highly skewed signals that influence the standard deviation.Finally, Table6 presents a brief recap of the steps necessary for reliableXApEn estimation.Table 6
Steps for reliableXApEn estimation.
Step
Explanation
(
1
)
Choose two time series and specify which one would be master series; defineN and m (preferably m=2); perform standard scoring of both series.
(
2
)
Find the differential time series both for master and for follower time series, and calculate their respective standard deviationsσDX and σDY.
(
3
)
CalculaterTH-A (see (18)).
(
4
)
CalculaterTH-X (see (17)).
(
5
)
Select reliability criterion; a strong criterion requires long data sets, so a weak criterion is suggested.
(
6
)
CalculaterXW (see (19)), weak criterion (preferably); or calculate rXS (see (20)), strong criterion.
(
7
)
EstimateXApEn(N, m, rXW).


## 3.1. Estimation Errors in Artificial Data

ForXApEn error estimation, master series X and a hundred of its followers Y were artificially generated for a range of (r,N,m,(τ=1)) parameters, and the same values exactly were calculated using the derived formulae. The results for all four data sets, due to the size, are presented within the Supplementary Material (available online at https://doi.org/10.1155/2017/8365685).Relative conditional probability errorsεim(r)/pim(r)[%] are plotted in Figure 3, for series X and Y both normally distributed and for τ=1. Satisfactory results are obtained for an extremely atypical threshold level r=1.50. For a typical threshold r=0.15, error is below 10% only for large N and m=1. The results are in accordance with reliability criteria N-m+1>(10or100)/pim(r) stating that N or pim(r) should be large. Indeed, pim(r) is defined as a product of m probabilities (14), so it is at its largest if m=1, while increased threshold r augments the number of template matches N1(i) (4), again enlarging pim(r).Figure 3
The relative conditional probability error[%] in artificial data for a typical threshold value r=0.15 (a), for a higher threshold r=0.60 (b), and for an unusually high threshold r=1.50 (c); the gray dashed horizontal line marks 10% error; results are presented in a log-log plane as mean±σ (standard deviation). Artificial master and follower series are both normally distributed.
(a)
(b)
(c)The reliability ofp^im(r) estimate is checked comparing the data length N-m+1 with the required minimum 10/pim(r). The percentage of test failures is presented in Figure 4, in an r-N plane, where 0% failures are presented by green color that gradually changes to yellow (100% failures). The white line marks classical r and N values, revealing that, for m=3 and 4, the estimates that implement classical parameters are unreliable. Nonlinear interrelation of r, N, and m in Figure 4 also points out that the choice of parameters that would yield a stable working point is not a straightforward task.Figure 4
Percentage of estimatesp^im(r) in artificial data that fail to satisfy the criterion N-m+1≥10/p^im(r)[%], presented in threshold-length (r-N) plane, for m=1,2,3,4; failures are presented by color changes from green (0% failures) to yellow (100% failures); (a) X and Y normally distributed; (b) X and Y uniformly distributed; white horizontal and vertical lines cut each plane at typical parameter values r=0.2 and N=1000.
(a)
(b)The estimation error of summandΦ^ (12) exposes two problems. The first one is a lack of matches (p^im(r)=0) that induce logarithms of zero in (6) which is, as shown in Figure 5, a subset of all unreliable probabilities from Figure 4. The second problem is an error accumulation (12).Figure 5
Number ofp^im(r)=0 estimates [%] in artificial data presented in threshold-length (r-N) plane for m=3 and 4; color changes from green (0% zero matches) to yellow (100% zero matches); (a) X and Y normally distributed; (b) X and Y uniformly distributed.
(a)
(b)Figure6 shows a line graph of summand Φ and its changes with respect to series length N if it is evaluated by a formula (exact values), if it is estimated in a traditional way, and if a correction (9) is applied. This correction does attenuate estimation errors, even for high m values, but it cannot be used alone, as it removes a subset of the problem only. The convergence of summand Φ^ towards its true value depends on data type: uniform data achieve the correct limit for shorter lengths N. The convergence is also m-dependent, requiring different data lengths to achieve the same limit. Since X^ApEn is defined as a difference of Φ^ values (7), the fact that each summand Φ^ adds its own and independent estimation error is a further source of X^ApEn instability.Figure 6
Summand Φ forr=0.3: estimated omitting the zero probability, estimated and corrected, and calculated exactly using a formula; (a) X and Y normally distributed; (b) X and Y uniformly distributed.
(a)
(b)Figure7 shows two-dimensionalXApEn profiles. Values calculated using the derived formulae ((a) and gray lines in (b) and (c)) verify the known fact thatXApEn linearly decreases with the logarithm of threshold value [5]. The major conclusion is that the trueXApEn values are influenced neither by time series length N nor by template length m. Frequent notions thatXApEn monotonously increases with N until a plateau is reached or that its value decreases with m are a consequence of probability underestimation that decreases with series length and eventually disappears. The panels in the second row of (a) show corrected X^ApEn estimates that are in good accordance with true values for m=1, but showing increased discrepancies for m=2 and 3. If no correction is applied (the third row of (a)),    X^ApEn substantially diverges from the correct values. (b) and (c) show X^ApEn as a line graph. The discrepancy between the true and estimated values points out that generally accepted guidelines for parameter choice do not offer a safe X^ApEn estimation.Figure 7
XApEn of artificial data with normal distribution in threshold-length (r-N) plane, for m=1,2, and 3. The three rows of (a) show, respectively, calculated (true)XApEn,XApEn estimated with correction, andXApEn estimated omitting log⁡(0). Dark green color corresponds toXApEn estimated within the unstableN-r-m parameter region. Light green indicates lowXApEn, while yellow color indicates highXApEn. Horizontal cuts at N=2000 are presented as line graph in (b); vertical cuts at r=0.3 are presented as line graph in (c).
(a)
(b)
(c)


## 3.2. Stable Working Point in Cardiovascular Environment

To verify the accuracy of the derived equation (17), we followed the recommendations from [15, 16], where their equations are validated with respect toApEn relative error. For each one of SBP-PI (or R-R) and PI (or R-R)-SBP signal pairs, of animals (or volunteers), a threshold rTH-X was evaluated according to (17) and used as a parameter forXApEn(rTH-X) estimation. The true maximalXApEn(rmax) value is provided by a time-consuming procedure of finding a threshold profile. The corresponding relative error εXA_max is presented in Table 2. The results estimated from signals recorded in human volunteers are presented in the italic part of Table 2. The relative error is below 2%, except for very short data lengths (N=300) and m=1 when it reaches 6%; unfortunately, this corresponds to the ambulatory records that cannot be too long.Table 2
Relative error ofXApEn(rTH-X) with respect to XApEn(rmax) in real signals [%].
N
εXA_max=XApEn(rTH-X)-XApEn(rmax)/XApEn(rmax)⋅100
X
=
P
I, Y=SBP
X
=
S
B
P, Y=PI
m = 1
m = 2
m = 3
m = 4
m = 1
m = 2
m = 3
m = 4
300
5.16 ± 0.40
1.25 ± 0.28
1.61 ± 0.21
5.51 ± 0.75
6.40 ± 0.61
2.42 ± 0.49
2.85 ± 0.65
5.60 ± 0.88
500
2.57 ± 0.25
1.13 ± 0.21
2.61 ± 0.51
2.91 ± 0.57
3.61 ± 0.39
1.74 ± 0.29
1.68 ± 0.34
2.83 ± 0.60
750
1.57 ± 0.18
0.90 ± 0.21
1.93 ± 0.49
2.06 ± 0.44
2.09 ± 0.26
1.09 ± 0.25
1.56 ± 0.28
1.62 ± 0.43
2000
0.64 ± 0.17
1.07 ± 0.26
1.50 ± 0.45
1.83 ± 0.54
0.49 ± 0.16
1.43 ± 0.29
1.96 ± 1.40
2.62 ± 0.47
2500
0.46 ± 0.07
0.88 ± 0.14
1.29 ± 0.29
1.61 ± 0.42
0.32 ± 0.05
1.01 ± 0.22
1.45 ± 0.29
2.99 ± 0.47
3000
0.80 ± 0.16
0.91 ± 0.18
0.99 ± 0.25
1.22 ± 0.26
0.33 ± 0.08
0.99 ± 0.18
0.94 ± 0.16
2.03 ± 0.38
3500
0.67 ± 0.12
1.04 ± 0.19
1.89 ± 0.67
2.68 ± 1.28
0.52 ± 0.14
1.41 ± 0.33
1.58 ± 0.47
2.64 ± 0.58
4000
1.03 ± 0.20
1.65 ± 0.28
1.58 ± 0.26
1.15 ± 0.23
1.08 ± 0.25
1.22 ± 0.23
1.09 ± 0.24
1.79 ± 0.45
4500
0.68 ± 0.18
1.41 ± 0.39
1.72 ± 0.59
2.80 ± 0.73
0.94 ± 0.18
1.20 ± 0.30
1.90 ± 0.92
2.28 ± 1.33
5000
1.59 ± 0.36
1.57 ± 0.44
1.73 ± 0.61
2.26 ± 0.75
0.42 ± 0.08
2.20 ± 0.39
1.74 ± 0.56
1.64 ± 0.40
Results are presented as mean ± SEM. Results of 41 human volunteers are shown in the first three italic rows. Results of 25 rats in different conditions (total of 50 signals) are presented in roman rows.Table 3
Relative error ofApEn(rTH-A) with respect to ApEn(rmax) [%] in real signals.
N
εA_max=ApEn(rTH-A)-ApEn(rmax)/ApEn(rmax)⋅100
X
=
P
I
X
=
S
B
P
r
TH-A [15]
r
TH-A [16]
r
TH-A [15]
r
TH-A [16]
2000
3.46 ± 0.53
4.07 ± 0.54
1.95 ± 0.33
2.33 ± 0.35
3000
3.23 ± 0.54
3.82 ± 0.56
2.30 ± 0.35
2.91 ± 0.48
4000
3.72 ± 0.65
4.40 ± 0.67
2.86 ± 0.48
3.33 ± 0.51
5000
3.75 ± 0.67
6.45 ± 1.26
2.68 ± 0.43
5.63 ± 1.29
Error is evaluated according to the guidelines from [15] and from [16]. Results are presented as mean ± SEM, averaged over 100 time series: 50 data series of Wistar and borderline hypertensive rats in baseline and in stress, repeated for m=2 and for m=3.For the sake of comparison, we estimatedApEn error from our data as well. We applied rTH-A calculated according to [15] (see (18)) and according to [16]. The rTH-A relative error is shown in Table 3, providing the following conclusions: (1) error level is lower if the equations from [15] are used and (2) comparison of results in Tables 2 and 3 verifies that our expressions (17) are equal in performance to expressions (18) derived by Lu et al. in [15, 16].The valuerTH-X cannot be applied, since the percentage of the probabilities that fail the weak and strong criteria is almost 100% (Figure 8). But this threshold reflects the nature of data and it can be used as a reference for thresholds rXW and rXS. Their implementation reduces the failures, as shown in Figure 8.Figure 8
Unreliable estimated probabilitiesp^im(r)[%] in real signals; (a) r∈{rXW,rTH-X} and weak criterion is applied; (b) r∈{rXS,rTH-X} and strong criterion is applied. The gray dashed horizontal line marks 10% of the unreliable probability estimates. The left part of each panel shows the results of Wistar and borderline hypertensive rats at baseline and stress; the right part of each panel shows the results of volunteers.
(a)
(b)The actual threshold values calculated from real signals are shown in Table4. For the sake of comparison, threshold values for normally distributed artificial data are also plotted (Figure 9). These thresholds show that the traditional guidelines (0.15, 0.2, and 0.35) considerably underestimate the actual requirements for the reliableXApEn estimation. Even if m=1, the threshold required for strong reliability may exceed r=3. Although it was speculated [1] that detailed system information would be lost if r values were larger than 0.1·σ or 0.2·σ (forApEn), meaning that if the threshold exceeds these values a template would always match its follower, in investigated cardiovascular signals, this assumption does not hold as the evaluated high threshold levels are considerably below the amplitude range of the observed data.Table 4
Threshold values obtained from real signals.
(a)
ThresholdrTH-X for which XApEn reaches its maximal value (17)
N
X
=
P
I, Y=SBP
X
=
S
B
P, Y=PI
m = 1
m = 2
m = 3
m = 4
m = 1
m = 2
m = 3
m = 4
300
0.05 ± 0.01
0.24 ± 0.01
0.44 ± 0.01
0.62 ± 0.01
0.05 ± 0.00
0.22 ± 0.01
0.40 ± 0.01
0.55 ± 0.02
500
0.05 ± 0.01
0.21 ± 0.00
0.39 ± 0.01
0.53 ± 0.00
0.04 ± 0.00
0.19 ± 0.00
0.35 ± 0.01
0.47 ± 0.01
750
0.04 ± 0.01
0.18 ± 0.04
0.35 ± 0.01
0.47 ± 0.01
0.04 ± 0.01
0.16 ± 0.00
0.31 ± 0.01
0.41 ± 0.01
2000
0.05 ± 0.01
0.19 ± 0.00
0.37 ± 0.01
0.51 ± 0.01
0.04 ± 0.00
0.14 ± 0.00
0.27 ± 0.00
0.60 ± 0.01
2500
0.05 ± 0.01
0.18 ± 0.00
0.35 ± 0.01
0.48 ± 0.01
0.04 ± 0.00
0.13 ± 0.00
0.25 ± 0.01
0.34 ± 0.01
3000
0.04 ± 0.01
0.18 ± 0.00
0.33 ± 0.01
0.46 ± 0.01
0.03 ± 0.00
0.13 ± 0.00
0.24 ± 0.00
0.33 ± 0.05
3500
0.04 ± 0.01
0.17 ± 0.00
0.32 ± 0.01
0.45 ± 0.01
0.03 ± 0.00
0.12 ± 0.00
0.23 ± 0.00
0.31 ± 0.05
4000
0.04 ± 0.01
0.17 ± 0.00
0.31 ± 0.01
0.44 ± 0.01
0.03 ± 0.00
0.17 ± 0.00
0.22 ± 0.00
0.30 ± 0.05
4500
0.04 ± 0.01
0.11 ± 0.00
0.30 ± 0.01
0.43 ± 0.01
0.03 ± 0.00
0.11 ± 0.00
0.21 ± 0.00
0.29 ± 0.01
5000
0.04 ± 0.01
0.16 ± 0.00
0.29 ± 0.01
0.42 ± 0.01
0.03 ± 0.00
0.11 ± 0.00
0.20 ± 0.00
0.29 ± 0.01
(b)
ThresholdrXW for which at least 95% of the probabilities satisfy N-m+1>10/p^im(rXW) (19)
N
X
=
P
I, Y=SBP
X
=
S
B
P, Y=PI
m = 1
m = 2
m = 3
m = 4
m = 1
m = 2
m = 3
m = 4
300
0.46 ± 0.00
0.84 ± 0.01
1.24 ± 0.01
1.62 ± 0.01
0.45 ± 0.00
0.82 ± 0.01
1.20 ± 0.01
1.55 ± 0.02
500
0.45 ± 0.00
0.81 ± 0.00
1.19 ± 0.01
1.53 ± 0.01
0.44 ± 0.00
0.78 ± 0.00
1.15 ± 0.01
1.47 ± 0.01
750
0.38 ± 0.00
0.69 ± 0.00
1.03 ± 0.01
1.32 ± 0.01
0.38 ± 0.00
0.67 ± 0.00
0.99 ± 0.01
1.26 ± 0.01
2000
0.39 ± 0.00
0.70 ± 0.00
1.05 ± 0.01
1.36 ± 0.01
0.37 ± 0.00
0.65 ± 0.00
0.95 ± 0.01
1.21 ± 0.01
2500
0.35 ± 0.00
0.63 ± 0.00
0.95 ± 0.01
1.23 ± 0.01
0.34 ± 0.00
0.58 ± 0.00
0.85 ± 0.01
1.09 ± 0.01
3000
0.34 ± 0.00
0.63 ± 0.00
0.93 ± 0.01
1.22 ± 0.01
0.33 ± 0.00
0.59 ± 0.00
0.84 ± 0.00
1.08 ± 0.01
3500
0.34 ± 0.00
0.62 ± 0.00
0.92 ± 0.01
1.20 ± 0.01
0.33 ± 0.00
0.57 ± 0.00
0.83 ± 0.00
1.06 ± 0.01
4000
0.34 ± 0.00
0.62 ± 0.00
0.91 ± 0.01
1.19 ± 0.01
0.33 ± 0.00
0.57 ± 0.00
0.82 ± 0.00
1.05 ± 0.01
4500
0.34 ± 0.00
0.61 ± 0.00
0.90 ± 0.01
1.18 ± 0.01
0.33 ± 0.00
0.56 ± 0.00
0.81 ± 0.00
1.04 ± 0.01
5000
0.34 ± 0.00
0.61 ± 0.00
0.89 ± 0.01
1.17 ± 0.01
0.33 ± 0.00
0.56 ± 0.00
0.80 ± 0.00
1.04 ± 0.01
(c)
ThresholdrXS for which at least 95% of the probabilities satisfy N-m+1>100/p^im(rXS) (20)
N
X
=
P
I, Y=SBP
X
=
S
B
P, Y=PI
m = 1
m = 2
m = 3
m = 4
m = 1
m = 2
m = 3
m = 4
300
1.42 ± 0.00
2.21 ± 0.02
3.28 ± 0.04
4.37 ± 0.09
1.40 ± 0.01
2.12 ± 0.02
3.07 ± 0.05
3.99 ± 0.10
500
1.08 ± 0.00
1.73 ± 0.14
2.60 ± 0.04
3.41 ± 0.01
1.06 ± 0.00
1.65 ± 0.01
2.41 ± 0.00
3.06 ± 0.07
750
0.96 ± 0.00
1.54 ± 0.01
2.28 ± 0.03
2.93 ± 0.05
0.96 ± 0.00
1.47 ± 0.01
2.13 ± 0.03
2.64 ± 0.04
2000
0.91 ± 0.00
1.50 ± 0.01
2.29 ± 0.03
3.06 ± 0.05
0.88 ± 0.00
1.33 ± 0.01
1.89 ± 0.02
2.29 ± 0.00
2500
0.90 ± 0.00
1.46 ± 0.01
2.21 ± 0.03
2.93 ± 0.05
0.88 ± 0.00
1.31 ± 0.01
1.82 ± 0.02
2.20 ± 0.03
3000
0.89 ± 0.00
1.44 ± 0.00
2.14 ± 0.03
2.83 ± 0.04
0.87 ± 0.00
1.29 ± 0.01
1.77 ± 0.02
2.14 ± 0.03
3500
0.89 ± 0.00
1.41 ± 0.01
2.08 ± 0.03
2.76 ± 0.04
0.87 ± 0.00
1.27 ± 0.01
1.72 ± 0.02
2.07 ± 0.02
4000
0.88 ± 0.00
1.40 ± 0.01
2.05 ± 0.02
2.72 ± 0.04
0.86 ± 0.00
1.25 ± 0.01
1.67 ± 0.02
2.00 ± 0.03
4500
0.88 ± 0.00
1.39 ± 0.01
2.00 ± 0.02
2.66 ± 0.04
0.86 ± 0.00
1.24 ± 0.01
1.63 ± 0.02
1.95 ± 0.02
5000
0.88 ± 0.00
1.37 ± 0.01
1.97 ± 0.02
2.62 ± 0.04
0.86 ± 0.00
1.23 ± 0.01
1.60 ± 0.00
1.93 ± 0.02
Results are presented as mean ± SEM, averaged over 50 time series of both Wistar and borderline hypertensive rats, both in baseline and in stress, and averaged over 41 time series of volunteers; results from human volunteers are in italics.Figure 9
Calculated threshold values for artificial data, with signalsX and Y both normally distributed. (a) Threshold for whichXApEn gets maximal value (see (17)). (b) Threshold for which at least 95% of the probabilities satisfy weak criterion (see (19)). (c) Threshold for which at least 95% of the probabilities satisfy strong criterion (see (20)). The gray horizontal dashed line corresponds to r=0.3.
(a)
(b)
(c)The developed sets of threshold expressions are primarily intended for cardiovascular signal analysis. The artificial signals used as well are uniform on segment (with bounded amplitudes) and Gaussian (with increased incidence of small amplitudes, but unbounded). Both artificial distributions are symmetric. The amplitudes of real signals (PI, SBP, and R-R) are bounded by physiological constraints. Signal asymmetry is expressed by skewness factor, shown in Table5. Cardiovascular data can be skewed to the left and to the right, so Table 5 presents the average of absolute skewness values. SBP and PI signals exhibit a level of asymmetry comparable to Rayleigh distribution (skewness ≈ 0.631), but below the skewness of exponential distribution (skewness ≈ 2). Threshold evaluation in signals with highly skewed distribution would be a subject for follow-up studies. In particular, it would be precious to explore an important class of interspike interval time series, with distributions that are close to exponential. A memoryless property of exponential distribution might affect the values of differential time series standard deviation σD, necessary for (17) and (18).Table 5
Level of asymmetry of real signals (skewness).
Rats
Volunteers
SBP
PI
SBP
R-R
Absolute skewness
0.64 ± 0.86
0.48 ± 1.02
0.42 ± 0.54
0.35 ± 0.33
Portion of signals with negative skewness[%]
36%
24%
46%
60%
Absolute skewness is expressed asmean±σ; each signal set comprises a couple of highly skewed signals that influence the standard deviation.Finally, Table6 presents a brief recap of the steps necessary for reliableXApEn estimation.Table 6
Steps for reliableXApEn estimation.
Step
Explanation
(
1
)
Choose two time series and specify which one would be master series; defineN and m (preferably m=2); perform standard scoring of both series.
(
2
)
Find the differential time series both for master and for follower time series, and calculate their respective standard deviationsσDX and σDY.
(
3
)
CalculaterTH-A (see (18)).
(
4
)
CalculaterTH-X (see (17)).
(
5
)
Select reliability criterion; a strong criterion requires long data sets, so a weak criterion is suggested.
(
6
)
CalculaterXW (see (19)), weak criterion (preferably); or calculate rXS (see (20)), strong criterion.
(
7
)
EstimateXApEn(N, m, rXW).


## 4. Conclusion

This paper performed an analysis ofXApEn parameter space and proposed methods for adjusting both the procedure and its parameter choice. Artificial environment with known data distribution allowed exact mathematical evaluation of referenceXApEn values and their comparison to the estimated ones. The comparison indicated vast regions inXApEn parameter space (threshold r, data length N, and template length m) with unstable estimation of bothXApEn and its constitutive components. It is also shown that parameter interrelationship is nonlinear and depends on data type. An established relationship between the theory of binary probability estimation [27] and template matching procedure proved that the source of instability is unreliable probability estimation, which is a consequence of working within the unstable region. StableXApEn estimates should remain constant for increased data lengths and the frequently observed monotonous increase with N is just another form of instability. It is also shown that the traditional guidelines work within the regions of instability. The proposed correction attenuated but did not eliminate the problem.The solution is to target the parameters within the stable region of the parameter space and to perform the parameter choice jointly, due to the nonlinear properties of their relationship. Real environment does not offer an advantage of exact mathematical error analysis, so a set of parameters for stable and reliable estimation is ensured deriving a novel set of thresholdr expressions. The expressions take into account length N, template m, statistical data properties, and reliability criterion (weak or strong, [27]). The consistency of estimates is validated for typical cases of short records (ambulatory patients) and long records (small laboratory animals). It is shown that classical guidelines underestimate the threshold levels and data lengths required for consistent results.This analysis can be employed for all the entropy measures that rely on probability estimation. The proposedXApEn improvements may contribute to its wider implementation in the domain of related physiological processes such as baroreflex studies. In order to facilitate and promote the widerXApEn deployment, we have included the detailed instructions within the Supplementary Material.Although the method is validated in cardiovascular data, it can be equivalently applied for finding correct parameters forXApEn analysis of processes that are not necessarily biomedical, provided that their distribution is not too skewed. Further investigation would consider the signals with highly asymmetric distribution, for example, interspike interval time series, but also the signals with mismatched distribution where zero probability is not an error, but inherent signal property.



---

*Source: 8365685-2017-09-17.xml*