AJP - Heart AJP: Heart and Circulatory Physiology
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
 QUICK SEARCH:   [advanced]


     


Am J Physiol Heart Circ Physiol 278: H2039-H2049, 2000;
0363-6135/00 $5.00
This Article
Right arrow Abstract Freely available
Right arrow Full Text (PDF)
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Right arrow Citation Map
Services
Right arrow Email this article to a friend
Right arrow Similar articles in this journal
Right arrow Similar articles in ISI Web of Science
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via ISI Web of Science (58)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Richman, J. S.
Right arrow Articles by Moorman, J. R.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Richman, J. S.
Right arrow Articles by Moorman, J. R.
Vol. 278, Issue 6, H2039-H2049, June 2000

Physiological time-series analysis using approximate entropy and sample entropy

Joshua S. Richman1,2 and J. Randall Moorman1

1 Cardiovascular Division, Department of Internal Medicine, and Department of Molecular Physiology and Biological Physics, and Cardiovascular Research Center, University of Virginia Health Sciences Center, Charlottesville, Virginia 22908; and 2 Medical Automation Systems, Charlottesville, Virginia 22903


    ABSTRACT
TOP
ABSTRACT
INTRODUCTION
THEORY
RESULTS AND DISCUSSION
REFERENCES

Entropy, as it relates to dynamical systems, is the rate of information production. Methods for estimation of the entropy of a system represented by a time series are not, however, well suited to analysis of the short and noisy data sets encountered in cardiovascular and other biological studies. Pincus introduced approximate entropy (ApEn), a set of measures of system complexity closely related to entropy, which is easily applied to clinical cardiovascular and other time series. ApEn statistics, however, lead to inconsistent results. We have developed a new and related complexity measure, sample entropy (SampEn), and have compared ApEn and SampEn by using them to analyze sets of random numbers with known probabilistic character. We have also evaluated cross-ApEn and cross-SampEn, which use cardiovascular data sets to measure the similarity of two distinct time series. SampEn agreed with theory much more closely than ApEn over a broad range of conditions. The improved accuracy of SampEn statistics should make them useful in the study of experimental clinical cardiovascular and other biological time series.

probability; nonlinear dynamics


    INTRODUCTION
TOP
ABSTRACT
INTRODUCTION
THEORY
RESULTS AND DISCUSSION
REFERENCES

NONLINEAR DYNAMICAL ANALYSIS is a powerful approach to understanding biological systems. The calculations, however, usually require very long data sets that can be difficult or impossible to obtain. Pincus (21, 22) devised the theory and method for a measure of regularity closely related to the Kolmogorov entropy, the rate of generation of new information, that can be applied to the typically short and noisy time series of clinical data. This family of statistics, named approximate entropy (ApEn), is rooted in the work of Grassberger and Procaccia (10) and Eckmann and Ruelle (5) and has been widely applied in clinical cardiovascular studies (3, 6-8, 12-19, 24, 26, 29, 32-34).

The method examines time series for similar epochs: more frequent and more similar epochs lead to lower values of ApEn. Informally, given N points, the family of statistics ApEn(mrN ) is approximately equal to the negative average natural logarithm of the conditional probability that two sequences that are similar for m points remain similar, that is, within a tolerance r, at the next point. Thus a low value of ApEn reflects a high degree of regularity. Importantly, the ApEn algorithm counts each sequence as matching itself, a practice carried over from the work of Eckmann and Ruelle (5) to avoid the occurrence of ln (0) in the calculations. This step has led to discussion of the bias of ApEn (22, 23, 27). In practice, we find that this bias causes ApEn to lack two important expected properties. First, ApEn is heavily dependent on the record length and is uniformly lower than expected for short records. Second, it lacks relative consistency. That is, if ApEn of one data set is higher than that of another, it should, but does not, remain higher for all conditions tested (22). This shortcoming is particularly important, because ApEn has been repeatedly recommended as a relative measure for comparing data sets (22-24).

To reduce this bias, we have developed and characterized a new family of statistics, sample entropy (SampEn), that does not count self-matches. SampEn is derived from approaches developed by Grassberger and co-workers (2, 9-11). SampEn(mrN ) is precisely the negative natural logarithm of the conditional probability that two sequences similar for m points remain similar at the next point, where self-matches are not included in calculating the probability. Thus a lower value of SampEn also indicates more self-similarity in the time series. In addition to eliminating self-matches, the SampEn algorithm is simpler than the ApEn algorithm, requiring approximately one-half as much time to calculate. SampEn is largely independent of record length and displays relative consistency under circumstances where ApEn does not.

Cross-ApEn is a recently introduced technique for analyzing two related time series to measure the degree of their asynchrony (20, 28). Cross-ApEn is very similar to ApEn in design and intent, differing only in that it compares sequences from one series with those of the second. Because it does not compare a series with itself, bias from self-matches does not arise. A potential problem, however, remains in the necessity for each template to generate a defined, nonzero probability. Thus each template must find at least one match for m + 1 points, or a probability must be assigned to it according to a "correction" strategy. We tested the effect of two extremes of correction strategies on cross-ApEn analysis. We find that cross-ApEn analysis lacks relative consistency, and conclusions about relative synchrony of pairs of time series depend on the unguided selection of analysis schemes. Cross-SampEn, on the other hand, is defined as long as one template finds a match, and we find that cross-SampEn remains relatively consistent for conditions where cross-ApEn does not.


    THEORY
TOP
ABSTRACT
INTRODUCTION
THEORY
RESULTS AND DISCUSSION
REFERENCES

ApEn reports on similarity in time series. We employ the terminology and notation of Grassberger and Procaccia (10), Eckmann and Ruelle (5), and Pincus (21) in describing techniques for estimating the Kolmogorov entropy of a process represented by a time series and the related statistics ApEn and SampEn. The parameters N, m, and r must be fixed for each calculation. N is the length of the time series, m is the length of sequences to be compared, and r is the tolerance for accepting matches. It is convenient to set the tolerance as r × SD, the standard deviation of the data set, allowing measurements on data sets with different amplitudes to be compared. Throughout this work, all time series have been normalized to have SD = 1.

We proceed as follows: For a time series of N points, {u( j): 1 <=  j <=  N } forms the N - m + 1 vectors xm(i) for {i|1 <=  i <=  N - m  + 1}, where xm(i) = {u(i k): 0 <=  k <=  m - 1} is the vector of m data points from u(i) to u(i + m - 1). The distance between two such vectors is defined to be d[x(i), xj)] = max {|u(i + k- uj + k)|:0 <=  k <=  m - 1}, the maximum difference of their corresponding scalar components. Let Bi be the number of vectors xm( j) within r of xm(i) and let Ai be the number of vectors xm + 1j ) within r of xm + 1(i). Define the function Cmi(r) = (Bi)/(N - m + 1). In calculating Cmi(r), the vector xm(i) is called the template, and an instance where a vector xm( j) is within r of it is called a template match. Cmi(r) is the probability that any vector xm( j) is within r of xm(i). The function Phi m(r) = (N - m + 1)-1 Sigma N - m + 1i = 1 lnt[Cmi(r)] is the average of the natural logarithms of the functions Cmi(r). Eckmann and Ruelle (5) suggest approximating the entropy of the underlying process as limr right-arrow  0 limm right-arrow  infinity limNright-arrow infinity  [Phi m(r- Phi m + 1(r)]. Because of these limits, this definition is not suited to the analysis of the finite and noisy time series derived from experiments.

Pincus (21) saw that the calculation of Phi m(r- Phi m + 1(r) for fixed parameters m, r, and N had intrinsic interest as a measure of regularity and complexity. He defines the related parameter ApEn(mr) = limN right-arrow  infinity [Phi m(r- Phi m + 1(r)], which for finite data sets is estimated by the statistic ApEn(mrN ) = Phi m(r- Phi m + 1(r). Algebraic manipulation reveals that ApEn(mr,N) = (N - m + 1)-1 Sigma N - m + 1i = 1 ln [Cmi(r)] - (N - m)-1 Sigma N - mi = 1 ln[Cm + 1i(r)]. When N is large, ApEn(mrN ) is approximately equal to (N - m)-1 Sigma N - mi = 1 [-ln (Ai/Bi)], the average over i of the negative natural logarithm of the conditional probability that d[xm + 1(i), xm + 1( j)] <=  r given that d[xm(i), xm( j)] <=  r. The difference between ApEn(mr, N ) and this average logarithmic probability is less than (N - m + 1)-1 ln (N - m + 1), which is <0.05 for N - m + 1 >=  90 and <0.02 for N - m + 1 >=  283. They differ because there are N - m + 1 templates of length m, but only N - m templates of length m + 1 (27). Thus, although the quantity CmN - m + 1 (r) is defined, Cm + 1N - m + 1 (r) is not, simply because the vector um + 1(N - m + 1) does not exist. For practical purposes, ApEn can be thought of as the negative natural logarithm of the probability that sequences that are close for m points remain close for an additional point. Because conditional probabilities lie between 0 and 1, the parameter ApEn(mr) is a positive number of infinite range. For finite N, however, the largest possible value of ApEn(mrN ) is -Phi m + 1(r<=  -ln lfloor (N - m)-1⌋, so ApEn(mrN<=  ln (N - m).

ApEn(m, r, N) is biased and suggests more similarity than is present. It is important to note that ApEn takes a template-wise approach to calculating this average logarithmic probability, first calculating a probability for each template. The ApEn algorithm thus requires that each template contribute a defined, nonzero probability. This constraint is overcome by allowing each template to match itself. Formally, because d[xm(i), xm(i)] = 0 <=  r, the ApEn algorithm counts each template as matching itself, a practice we will refer to as self-matching. This ensures that the functions Cmi(r) are > 0 for all i, thereby avoiding the occurrence of ln (0) in the calculation. Thus ApEn(mrN ) is certain to be defined under all circumstances. Pincus and Goldberger (23, 27) note that, as a consequence of this practice, ApEn(mrN ) is a biased statistic. Formally, a statistic is biased if its expected value is not equal to the parameter it estimates. ApEn statistics are biased, because the expected value of ApEn(mrN ) is less than the parameter ApEn(mr) (27).

To discuss the bias caused by including self-matches, let us redefine the conditional probability associated with the template xi(m) by letting Bi denote the number of vectors xm( j) with j not equal  i, such that d[xm + 1(i), xm + 1( j)] <=  r by Ai. The ApEn algorithm thus assigns to the template xi(m) a biased conditional probability of (Ai + 1)/(Bi + 1), which is always greater than the unbiased Ai/Bi. In the limit as N approaches infinity, Ai and Bi will generally be large, making the biased and unbiased probabilities asymptotically equivalent. Therefore, this bias is evident only for the analysis of finite data sets and is a characteristic of the statistic ApEn(mrN ), rather than the parameter ApEn(m, r). For a finite N, however, the result is that ApEn(mrN ) is biased toward lower values of ApEn and returns values below those predicted by theory.

The largest deviation occurs when a large proportion of templates have Bi = Ai = 0, since these templates are assigned a conditional probability of 1, corresponding to perfect order. Furthermore, the difference between the biased and unbiased conditional probabilities assigned to individual templates makes the calculation sensitive to record length in a way that depends on the conditional probability. Suppose that the unbiased conditional probability is known and denote it by CP. For a given um(i) let Bi denote the number of template matches without counting self-matches. The original algorithm estimates CP as (1 + Ai)/(1 + Bi) = (1 + Bi × CP)/(1 + Bi). The fractional error of this relative to CP is Err = {[(1 + Bi × CP)/ (1 + Bi)] - CP} /CP). To find the value of Bi necessary to keep the fractional error below a threshold Errmax, we isolate Bi from the inequality Err <=  Errmax, yielding Bi >=  [1 - CP(Errmax + 1)]/(CP × Errmax). For independent, identically distributed (iid) random numbers, obtaining Bi matches of length m requires a data set containing, on average, Bi/(CP)m templates. For iid random numbers and m = 2, estimating CP = 0.368 [ApEn(mrN ) = 1] within Errmax = 0.05 requires Bi >=  33 and a data set of >240 points. Estimating CP = 0.135 [ApEn(m, rN ) = 2] with similar resolution requires Bi >=  127 and >6,900 points.

The most straightforward way to eliminate the bias would be to remove self-matching from the ApEn algorithm, leaving it otherwise unaltered. However, without the inclusion of self-matches, the ApEn algorithm is not defined unless Cm + 1i(r) > 0 for every i. Removing self-matches would make ApEn statistics highly sensitive to outliers; if there were a single template that matched no other vector, ApEn could not be calculated because of the occurrence of ln (0). For the set of uniform random numbers shown in Fig. 1A, without self-matches, ApEn(1, rN ) would not have been defined for r <=  0.63 and N < 4,000; ApEn(2, rN ) would have required r >=  1 for N < 4,000. Thus, for many practical applications, self-matches cannot simply be excluded when ApEn statistics are calculated.


View larger version (18K):
[in this window]
[in a new window]
 
Fig. 1.   Sample entropy (SampEn) and approximate entropy (ApEn) of random numbers with a uniform distribution. A: frequency histogram of 1,000 numbers. Inset: 25-point excerpt from time series. B: SampEn and ApEn as functions of r(m = 2). Straight line, calculated theoretical values for ApEn and SampEn parameters. In this special case of uniformly distributed random numbers, their theoretical values coincide. Values of r are displayed logarithmically, so that theoretical values appear linear. C: SampEn and ApEn as functions of N(m = 2, r = 0.2). In B and C, confidence intervals for SampEn statistics are displayed as error bars. Straight line is theoretically predicted value of parameter SampEn(2,0.2). Sets of uniformly distributed random numbers were generated using a minimal standard number generator with added random shuffling (30), and all sets passed a runs test for random arrangement (1). For data shown above and for similar tests in which random data with different probabilistic distributions were used, theoretical values were calculated from expressions for SampEn(mr) and ApEn(mr) by use of trapezoidal numerical integration. We verified our calculation of ApEn statistics by comparing our results with published values for Henon and logistic map data (21).

Can the bias be corrected? It is suggested that this bias can be reduced with a family of estimators epsilon  by defining Cmi,epsilon (r) = (N - m + 1)-1 {epsilon  + number of j not equal  i such that d[xm(i), xm( j)] <=  r} (27). Then Phi mepsilon (r) = (N - m + 1)-1 Sigma N - m + 1i = 1 Cmi,epsilon (r) and ApEnepsilon (m, rN ) = Phi mepsilon (r- Phi m + 1epsilon (r). For large N, this is very close to estimating the conditional probabilities by (epsilon  + Ai)/(epsilon  + Bi). When epsilon  < 1, the error due to counting self-matches is reduced, but this algorithm would still report a conditional probability of 1 in the event that only self-matches are encountered. Another strategy would be to define ApEnepsilon 2,epsilon 1(mrN ) = Phi mepsilon 2(r) - Phi m + 1epsilon 1(r) where epsilon 1 and epsilon 2 are small and epsilon 1 <=  epsilon 2. This is similar to estimating the conditional probabilities by (epsilon 1 + Ai)/ (epsilon 2 + Bi). Thus a template reporting only self-matches would be assigned a probability of epsilon 1/epsilon 2. This remains a biased estimate of ApEn, with bias toward -ln (epsilon 1/epsilon 2) when a large proportion of templates report only self-matches. If epsilon 1/epsilon 2 = (N - m + 1)-1 and Ai = Bi = 0, the bias is toward the lowest observable probability. This is the opposite bias of the original formulation of ApEn, where the bias is toward the highest possible probability. Still another approach would be to use the estimators epsilon 1 and epsilon 2 but incorporate them into the calculation only when Ai or Bi = 0, respectively. This would estimate the probabilities as Ai/Bi and would simply redefine Ai = epsilon 1 when Ai = 0 and Bi epsilon 2 when Bi = 0.

These approaches reduce the bias in the estimation of the individual conditional probabilities and ensure that perfect order would not be reported where none had been detected. None of the corrections eliminates bias; the bias is minimized only if epsilon 1/epsilon 2 = CP, but CP is not known beforehand. Thus no family of estimators epsilon  that minimizes bias can be chosen a priori.

SampEn statistics have reduced bias. We developed SampEn statistics to be free of the bias caused by self-matching. The name refers to the applicability to time series data sampled from a continuous process. In addition, the algorithm suggests ways to employ sample statistics to evaluate the results, as explained below.

There are two major differences between SampEn and ApEn statistics. First, SampEn does not count self-matches. We justified discounting self-matches on the grounds that entropy is conceived as a measure of the rate of information production (5), and in this context comparing data with themselves is meaningless. Furthermore, self-matches are explicitly dismissed in the later work of Grassberger and co-workers (2, 9, 11). Second, SampEn does not use a template-wise approach when estimating conditional probabilities. To be defined, SampEn requires only that one template find a match of length m + 1.

We began from the work of Grassberger and Procaccia (10), who defined Cm(r) = (N - m + 1)-1 Sigma N - m + 1i = 1 Cmi(r), the average of the Cmi(r) defined above. This differs from Phi m(r) only in that Phi m(r) is the average of the natural logarithms of the Cmi(r). They suggest approximating the Kolmogorov entropy of a process represented by a time series by limr right-arrow  0 limn right-arrow  infinity limN right-arrow  infinity  - ln [Cm + 1(r)/ Cm(r)]; self-matches are counted, and Cm + 1(r)/Cm(r) = (N - m + 1) Sigma N - mi = 1 Ai/(N - m) Sigma N - m + 1i = 1 Bi.

In this form, however, the limits render it unsuitable for the analysis of finite time series with noise. We therefore made two alterations to adapt it to this purpose. First, we followed their later practice in calculating correlation integrals (2, 9, 11) and did not consider self-matches when computing Cm(r). Second, we considered only the first N - m vectors of length m, ensuring that, for 1 <=  i <=  N - m, xm(i) and xm + 1(i) were defined.

We defined Bmi(r) as (N - m - 1)-1 times the number of vectors xm( j) within r of xm(i), where j ranges from 1 to N - m, and j not equal  i to exclude self-matches. We then defined Bm(r) = (N - m)-1 Sigma N - mi = 1 Bmi(r). Similarly, we defined Ami(r) as (N - m - 1)-1 times the number of vectors xm + 1( j) within r of xm + 1(i), where j ranges from 1 to N - m ( j not equal  i), and set Am(r) = (N - m)-1 Sigma N - mi = 1 Ami(r). Bm(r) is then the probability that two sequences will match for m points, whereas Am(r) is the probability that two sequences will match for m + 1 points. We then defined the parameter SampEn(mr) = limN right-arrow  infinity {-ln [Am(r)/ Bm(r)]}, which is estimated by the statistic SampEn(mrN ) = -ln [Am(r)/Bm(r)]. Where there is no confusion about the parameter r and the length m of the template vector, we set B = {[(N - m - 1)(N - m)]/2} Bm(r) and A = {[(N - m - 1)(N - m)]/2} Am(r), so that B is the total number of template matches of length m and A is the total number of forward matches of length m + 1. We note that A/B = [Am(r)]/Bm(r)], so SampEn(mrN ) can be expressed as -ln (A/B).

The quantity A/B is precisely the conditional probability that two sequences within a tolerance r for m points remain within r of each other at the next point. In contrast to ApEn(mrN ), which calculates probabilities in a template-wise fashion, SampEn(mrN ) calculates the negative logarithm of a probability associated with the time series as a whole. SampEn(mrN ) will be defined except when B = 0, in which case no regularity has been detected, or when A = 0, which corresponds to a conditional probability of 0 and an infinite value of SampEn(mrN ). The lowest nonzero conditional probability that this algorithm can report is 2[(N - m - 1)(N - m)]-1. Thus, the statistic SampEn(mrN ) has ln (N - m) + ln (N - m - 1) - ln (2) as an upper bound, nearly doubling ln (N - m), the dynamic range of ApEn(mrN ).

Confidence intervals inform the implementation of SampEn statistics. SampEn is not defined unless template and forward matches occur and is not necessarily reliable for small numbers of matches. We have reviewed SampEn(mrN ) calculation as a process of sampling information about regularity in the time series and used sample statistics to inform us of the reliability of the calculated result. For example, say that we find B template matches, allowing for no more than B forward matches, A of which actually occur. We assign a value of 1 to the A forward matches and a value of 0 to the (B - A) potential forward matches that do not occur and compute the conditional probability measured by SampEn as the average of this sample of 0s and 1s. For operational purposes, we will assume that the sample averages follow a Student's td distribution, where d is the number of degrees of freedom. We can then say with 95% confidence that the "true" average conditional probability of the process is within SDt(B - 1,0.975)/<RAD><RCD><IT>B</IT></RCD></RAD> of the sample average, where SD is the sample standard deviation and tB - 1,0.975 is the upper 2.5th percentile of a t distribution with B - 1 degrees of freedom (31). The size of the confidence intervals depends on the number B and the number of forward matches. Informally, large confidence intervals around SampEn(mrN ) indicate that there are insufficient data to estimate the conditional probability with confidence for that choice of m and r. In addition, confidence intervals allow standard statistical tests of the significance of differences between data sets.

The confidence intervals for SampEn are displayed as error bars in the figures. For some small values of N and r, no value of SampEn is given. This indicates that B = 0, A = 0, or the confidence intervals extended to a probability of >1 or <0. In these cases, no value of SampEn(mrN ) can be assigned with confidence.

Cross-ApEn and cross-SampEn measure asynchrony. Cross-ApEn is a recently introduced technique for comparing two different time series to assess their degree of asynchrony or dissimilarity (20, 28). The definition of cross-ApEn is very similar to ApEn. Given two time series of N points {u( j): <=  j <=  N } and {v( j): 1 <=  j <=  N }, form the vectors xm(i) = {u(i + k): 0 <=  k <=  m - 1} and ym(i) = {v(i + k): 0 <=  k <=  m - 1}. The distance between two such vectors is defined as d[xm(i), ym( j)] = max {|u(i k- v( j + k)|: 0 <=  k <=  m - 1}. Define Cmi(r)(v∥u) as the number of ym( j) within r of xm(i) divided by (N - m + 1), then define Phi m(r)(v∥u) = (N - m + 1)-1 Sigma N - m + 1i = 1 ln [Cmi(r) (v∥u)], and cross-ApEn(m, rN )(v∥u) = Phi m(r)(v∥u- Phi m + 1(r)(v∥u). This is identical to the definition of the statistic ApEn, except the templates are chosen from the series u and compared with vectors from v. There is thus a directionality to this analysis, and we will call the series that contributes the templates the template series and the series with which they are compared the target series. We can then refer to cross-ApEn(mrN ) (target∥template).

We made two observations. First, because no template is compared with itself, there are no self-matches. Consequently, Cmi(r)(v∥u) can equal 0, and there is no assurance that cross-ApEn(mrN )(v∥u) will be defined. Second, there is a "direction dependence" of cross-ApEn analysis. We defined cross-SampEn to avoid both potential problems.

Cross-ApEn is not always defined. As noted above, cross-ApEn does not include self-matching and, thus, does not inherently suffer from the same bias as ApEn. A potential problem, however, remains in the necessity for each template to generate a defined, nonzero conditional probability. Thus each template must find at least one match for m + 1 points, or a probability must be assigned to it. No guidelines have been suggested for handling this potential difficulty. Cross-SampEn, on the other hand, requires only that one pair of vectors in the two series match for m + 1 points.

The family of MIX(P) stochastic processes (21) provided a testing ground for cross-ApEn. Informally, the MIX(P) time series of N points, where P is between 0 and 1, is a sine wave, where N × P randomly chosen points have been replaced with random noise. We calculated cross-ApEn(1, r, 250) for the pair [MIX(Q)∥MIX(P)] and its direction conjugate [MIX(P)∥MIX(Q)] for 16 realizations of each of the 6 combinations of P = 0.1, 0.2, 0.3 and Q = 0.5, 0.7 over a range of values of r from 0.01 to 1.0. Cross-ApEn(1, r, 250) [MIX(Q)∥MIX(P)] was not defined for any of the 96 pairs for r <=  0.16 and was defined for all of them only for r >=  0.50. Cross-ApEn (1, r, 250) [MIX(P)∥MIX(Q)] was not defined for any values of r <=  0.32 and was defined for all pairs only for r = 1.0.

To broaden the conditions for which cross-ApEn was defined, we introduced a correction factor into its algorithm. To avoid ln (0) whenever Cmi(r)(v∥u) = 0 or Cm + 1i(r)(v∥u) = 0, we redefined them to be positive and nonzero. Rather than include these factors into the calculation of each Cmi(r)(v∥u) and Cm + 1i(r)(v∥u), we minimized their impact on the overall calculation by including them only when needed to ensure that cross-ApEn is defined. The cost of this adjustment, however, is the introduction of bias. For this reason, we tested cross-ApEn with two different correction strategies. The first strategy, which we called bias 0, was very similar to self-matching; we simply set any Cmi(r)(v∥u) or Cm + 1i(r)(v∥u) = 1 if it would otherwise have been 0. Thus, if a template matched no other at all, it would be assigned a conditional probability of 1, as with the original description of ApEn. If, however, Cmi(r)(v∥u) not equal  0 but Cm + 1i(r)(v∥u) = 0, we redefined Cm + 1i(r)(v∥u) = (N - m)-1, so that the probability assigned would be the lowest observable, nonzero probability given the nonzero value of Cmi(r)(v∥u).

The second approach, which we called bias max, also only modified the functions Cmi(r)(v∥u) and Cm + 1i(r) (v∥u) that would otherwise have been 0. Here, Cmi(r) (v∥u) was redefined to be 1 when it would otherwise have been 0 and Cm + 1i(r)(v∥u) was redefined to be (N - m + 1)-1, as for bias 0.

The only difference between the two strategies is that bias max assigns to a template yielding no matches at all a probability of (N - m)-1, the lowest nonzero probability allowed by the length of the time series. Thus bias 0 sets the bias toward a cross-ApEn value of 0 in the absence of any matches, whereas bias max sets the bias toward the highest observable value of cross-ApEn.

Cross-ApEn is direction dependent; cross-SampEn is not. Because of the logarithms inside the summation, Phi m(r)(v∥u) will not generally be equal to Phi m(r)(u∥v). Thus cross-ApEn(mrN )(v∥u) and its direction conjugate cross-ApEn(mrN )(u∥v) are unequal in most cases.

In defining cross-SampEn, we set Bmi(r)(v∥u) as (N - m)-1 times the number of vectors ym( j) within r of xm(i), where j ranges from 1 to N - m. We then defined Bm(r)(v∥u) = (N - m)-1 Sigma N - mi = 1 Bmi(r)(v∥u). Similarly, we set Ami(r)(v∥u) as (N - m)-1 times the number of vectors ym + 1( j) within r of xm + 1(i), where j ranges from 1 to N - m. We then defined Am(r)(v∥u) = (N - m)-1 Sigma N - mi = 1 Ami(r)(v∥u). Finally, we set cross-SampEn(mrN )(v∥u) = -ln {[Am(r)(v∥u)]/ [Bm(r)(v∥u)]}. Examining this definition for direction dependence, we see that (N - m) Bmi(r)(v∥u) is the number of vectors from v within r of the ith template of the series u. Summing over the templates, we see that Sigma N - mi = 1 (N - m) Bmi(r)(v∥u) simply counts the number of pairs of vectors from the two series that match within r. The number of pairs that match is clearly independent of which series is the template and which is the target. Because the last summation is equal to (N - m)2 Bm(r)(v∥u), it follows that Bm(r)(v∥u) is also direction independent, implying that cross-SampEn(m, rN )(v∥u) = cross-SampEn(mrN )(u∥v). We further noted that cross-SampEn will be defined provided that Am(r)(v∥unot equal  0. Cross-SampEn, on the other hand, requires only that one pair of vectors in the two series match for m + 1 points.

We calculated cross-SampEn(1, r, 250) for the same realizations of the [MIX(Q)∥MIX(P)] pairs used above to test cross-ApEn and over the same range of r. In contrast to cross-ApEn(1, r, 250) [MIX(P)∥MIX(Q)] and cross-ApEn(1, r, 250) [MIX(Q)∥MIX(P)], cross-SampEn(1, r, 250) [MIX(P)∥MIX(Q)], which is identical to cross-SampEn(1, r, 250) [MIX(Q)∥MIX(P)], was defined for all 96 pairs over the entire range of r considered.

ApEn and SampEn can be calculated analytically for series of random numbers. ApEn and SampEn derive from formulas suggested to estimate the Kolmogorov entropy of a process represented by a time series. At their root, each is a measurement of the conditional probability that two vectors that are close to each other for m points will remain close at the next point. There are several models, including sets of iid random numbers, for which the theoretical values of the parameters ApEn(mr)(21) and SampEn(mr) can be calculated. We show here the case of uniform random numbers, for which the theoretical values of ApEn and SampEn are nearly identical.

The expected value of the key probability can be calculated analytically for series of iid numbers based only on their probabilistic distribution. The numbers' independence implies that the probability that two randomly selected sequences within rSD of each other for the first m points will remain within r at their next points is simply the probability that any two points will be within a distance rSD of each other. For random numbers with density function p(x) and standard deviation SD, the expression for this probability is int +infinity -infinity p(t) [int t + rSDt - rSD p(x) dx] dt. Because the parameter ApEn(mr) measures the average of the negative logarithm of this probability, it is expected to return a value of - int +infinity -infinity p(t) ln [int t + rSDt - rSD p(x) dx] dt (21). SampEn(mr), on the other hand, measures the logarithm of the conditional probability and is thus expected to return a value of - ln {int +infinity -infinity p(t) [int t + rSDt - rSD p(x) dx] dt}.

Because these expressions for the expected values of the parameters ApEn and SampEn depend solely on r and the probabilistic character of the data, uniform random numbers provide a benchmark for testing the estimation of ApEn over a range of parameters. In particular, ApEn and SampEn are expected to give identical results for uniformly distributed random numbers. Figure 1A shows a histogram of the random numbers used and an excerpt from the sequence (inset). Taking the derivative of -int +infinity -infinity p(t) ln [int t + rSDt - rSD p(x) dx] dt with respect to ln (r) reveals that, for small r (in this case r <=  1), ApEn(mr) decreases in proportion to ln (r) for small r. This result generalizes to random numbers with other distributions provided that their density functions are smoothly continuous and bounded. The expected values are shown as the straight line in Fig. 1B. For random numbers with other distributions, the theoretical values of ApEn and SampEn will differ slightly. This is explained by noting that because -ln is a convex function, Jensen's inequality (4) applied to each of the integral expressions implies that the ApEn is expected to return a value greater than or equal to SampEn.


    RESULTS AND DISCUSSION
TOP
ABSTRACT
INTRODUCTION
THEORY
RESULTS AND DISCUSSION
REFERENCES

SampEn agrees with theory more closely than ApEn. For most processes, ApEn statistics are expected to have two properties. First, the conditional probability that sequences within r of each other remain within r should decrease as r decreases and the criterion for matching becomes more stringent. In other words, ApEn(mrN ) should increase as r decreases (22, 27). This expected property is demonstrated in Fig. 1B by the straight line, which plots the theoretically predicted values of ApEn and SampEn. Second, ApEn should be independent of record length. The plot of theoretical values of ApEn and SampEn in Fig. 1B illustrates this expectation. Because of their similarity, we expect SampEn statistics to exhibit similar properties. It has been suggested that record lengths of 10m-20m should be sufficient to estimate Ap- En(mr) (27).

We tested ApEn and SampEn statistics on uniform, iid random numbers, because the results could be compared with the analytically calculated expected values. Figure 1, B and C, shows the performance of ApEn(2, rN ) and SampEn(2, rN ) on uniform iid random numbers. SampEn(2, rN ) very closely matches the expected results for r >=  0.03 and N >=  100 (Fig. 1, B and C), whereas ApEn(2, rN ) differs markedly from expectations for N < 1,000 and r < 0.2. Figure 2 shows SampEn and ApEn as functions of r (m = 2) for three sets of uniform random numbers consisting of 100, 5,000, and 20,000 points. SampEn statistics for r = 0.2 are in agreement with theory for much shorter data sets. We investigated the general applicability of SampEn and ApEn statistics to random numbers with other distributions. The analysis of numbers with Gaussian, exponential, and gamma -distributions with parameter lambda  = 1,2, ... , 10 gave results essentially identical to those shown in Fig. 1.


View larger version (12K):
[in this window]
[in a new window]
 
Fig. 2.   Effect of record length on SampEn and ApEn as a function of r with m = 2 for uniform random numbers. Confidence intervals of SampEn calculations are displayed as error bars. Plots are SampEn and ApEn as functions of r for sets of 100 (A), 5,000 (B), and 20,000 (C) points. Straight lines, theoretically predicted values of ApEn and SampEn statistics.

SampEn shows relative consistency where ApEn does not. A critically important expected feature of ApEn is relative consistency (22, 27). That is, it is expected that, for most processes, if ApEn(m1r1)(S<=  ApEn(m1r1)(T ), then ApEn(m2, r2)(S), <=  ApEn(m2r2)(T ). That is, if record S exhibits more regularity than record T for one pair of parameters m and r, it is expected to do so for all other pairs. Graphically, plots of ApEn as a function of r for different data sets should not cross over one another. The determination that one set of data exhibits greater regularity than another can be made only when this condition is met. We tested this expectation using 1,000-point realizations of the MIX(P) process, where the degree of order could be specified. Figure 3A shows that the ApEn statistics of the less-ordered MIX(0.9), which has, on average, few matches of length m for a given template (small Bi) for small r, rises as a function of r and crosses over a plot of ApEn statistics of the more-ordered MIX(0.1). For r < 0.05, one would conclude incorrectly that MIX(0.9) was more ordered than MIX(0.1). Thus relative consistency does not hold for ApEn statistics.


View larger version (14K):
[in this window]
[in a new window]
 
Fig. 3.   Improved relative consistency but residual bias in SampEn. A and B: relative consistency is maintained by SampEn but not by ApEn for MIX(0.1) and MIX(0.9) processes. A: ApEn statistics as a function of r(m = 2, N = 1000) for MIX(0.1) and MIX(0.9). Inset: 25 sample points of MIX(0.1) (top) and MIX(0.9) (bottom). B: SampEn statistics for MIX(0.1) and MIX(0.9). C: residual bias of SampEn(2,0.2,N ) statistics for short record lengths of independent identically distributed Gaussian numbers. For N = 4-102, SampEn(2,0.2,N ) was calculated for >= 105 sets of random numbers and averaged.

We investigated the mechanism responsible for this lack of relative consistency. Note that, for r = 0.5, ApEn statistics correctly distinguished MIX(0.1) and MIX(0.9). For this value of r, the MIX(0.1) data yielded an average of >46 matches per template and ~28 forward matches per template, whereas the MIX(0.9) data yielded ~37 and 10, respectively. These large numbers of matches render the bias insignificant; that is, the unbiased Ai/Bi (28/46 and 10/37) is not very different from the biased (Ai + 1)/(Bi + 1) (29/47 and 11/38). As r is decreased, however, the number of template matches decreased more for the MIX(0.9) data than for the MIX(0.1) data. This made the bias significant for MIX(0.9) data under conditions for which it was insignificant in the MIX(0.1) data. For example, when r = 0.05, ApEn(2, r, 1,000) of MIX(0.1) was 0.463, whereas the value for MIX(0.9) was 0.505. The spurious similarity of the ApEn statistics is due to bias. The MIX(0.1) data yielded ~27 matches per template and 22 forward matches, whereas the MIX(0.9) data yielded only 0.45 and 0.01, respectively. Thus more than one-half of the templates of the MIX(0.9) data matched no other templates and were assigned a conditional probability of 1. In this example, ApEn statistics lack relative consistency, because less-ordered data sets have fewer matching templates and are more vulnerable to the bias generated by self-matches. SampEn analysis of the same data, on the other hand, reports correctly over the whole range of r (Fig. 3B).

Are SampEn statistics relatively consistent? We have shown that SampEn statistics appear to be relatively consistent over the family of MIX(P) processes, whereas ApEn statistics are not. Although we believe that relative consistency should be preserved for processes for which probabilistic character is understood, we see no general reason why ApEn or SampEn statistics should remain relatively consistent for all time series and all choices of parameters.

We propose a general, but by no means exhaustive, explanation for this phenomenon. SampEn is, in essence, an event-counting statistic, where the events are instances of vectors being similar to one another. When these events are sparse, the statistics are expected to be unstable, which might lead to a lack of relative consistency. Recall that SampEn(mrN ) is less than or equal to ln (B), the natural logarithm of the number of template matches. Suppose SampEn(mrN )(S) < SampEn(mrN ) (T ) and that the number of T's template matches, BT, is less than the number of S's template matches, BS, which would be consistent with T displaying less order than S. Provided that AT and AS, the number of forward matches, are relatively large, both SampEn statistics will be considerably lower than their upper bounds. As r decreases, BT and AT are expected to decrease more rapidly than BS and AS. Thus, as BT becomes very small, SampEn(mrN )(T ) will begin to decrease, approaching the value ln (BT), and could cross over a graph of SampEn(mrN )(S), where or while BS is still relatively large. Furthermore, as the number of template matches decreases, small changes in the number of forward matches can have a large effect on the observed conditional probability. Thus the discrete nature of the SampEn probability estimation could lead to small degrees of crossover and intermittent failure of relative consistency, and we cannot say that SampEn will always be relatively consistent. We have shown, however, that SampEn is relatively consistent for conditions where ApEn is not, and we have not observed any circumstance where ApEn maintains relative consistency and SampEn does not.

One source of the residual bias in SampEn is correlation of templates. Although more consonant with theory than ApEn, we found that SampEn statistics deviated from predictions for very short data sets. For 105 sets of Gaussian random numbers with m = 2, and r = 0.2, we found that the deviation was <3% for record lengths >100 but as high as 35% for sets of 15 points. Figure 3C shows the biased results of SampEn(2, 0.2, N ) for the range of 4 <=  N <=  100. We suspected that the integral expressions for the parameters ApEn(mr) and SampEn(mr) could not be used as expected values of the statistics ApEn(mrN ) and SampEn(mrN ) under all conditions, because the expressions relied on the assumption that the templates were independent of one another. As N decreases and m increases, however, a larger proportion of templates are comprised of overlapping segments of the record and are thus not independent. Because of this correlation, results might deviate from these predictions for short data sets. We thus tested the hypothesis that the majority of the bias results from nonindependence of the templates.

One way to test the hypothesis is to compare observed values of SampEn with those obtained from a model accounting for template correlation. The hypothesis predicts that the values should match. We tested the simplest case of m = 2 and N = 4, where a data set can be represented by {wxyz}, so that there are exactly two vectors of length m = 3, (wxy) and (xyz). If (wx) is close to (xy), it stands to reason that (wxy) will have a higher than expected probability of being close to (xyz). Formally, the conditional probability that (wxy) and (xy, z) will be close given (wx) and (xy) are close is [int infinity -infinity p(w) (int w + rSDw - rSD p(x) {int x + rSDx - rSD p( y) [int y + rSDy - rSD p(z) dz] dy} dx) dw]/ (int infinity -infinity p(w) {int w + rSDw - rSD p(x) [int x + rSDx - rSD p( y) dy] dx} dw), where p(z) is the Gaussian probability density function. For r = 0.2, this was evaluated by numerical integration and found to be 0.137. For 106 sets of 4 iid Gaussian numbers, we calculated SampEn(2, 0.2, 4) to be 0.140. Thus the observed and expected results nearly agree, indicating that bias due to nonindependence of templates accounts for most of the observed deviation from the conditional probability of 0.112 expected for independent templates in this case.

A second way to test the hypothesis is to calculate SampEn statistics for one set of data under two conditions: with and without overlapping templates. The hypothesis predicts that the results should not match. For this test, we chose the case of m = 2 and N = 6, the shortest record containing two nonoverlapping templates of length m + 1 = 3. We can represent each set as {abcx, yz}. For 106 sets of six Gaussian random numbers, we calculated the conditional probability that (abc) was within r of (xyz) given that their first two points were close, thus calculating the probability for pairs of disjoint templates. The result was 0.111, very close to the expected 0.112 for independent templates. For the same number sets, we then calculated the average value of SampEn(2, 0.2, 6), and the result was 0.094. Thus the two results do not match, in support of the hypothesis. We conclude from this analysis that the statistics SampEn(mrN ) are not completely unbiased under all conditions and that the bias of SampEn for very small data sets is largely due to nonindependence of templates.

One method for removing this bias would be to partition the time series {uj|1 <=  j <=   N } into the m + 1 sets of neighboring, disjoint vectors of length m + 1 Xi = {[ui + k(m + 1), u i + k(m + 1) + 1, ... , u i + k (m + 1) + m]: <=  k <=  [N - (m + i)]/(m + 1)}, where i, the initial point of the first template, ranges from 1 to m + 1. The conditional probability that vectors close for m points remain close at the next point would be calculated for each of the sets of vectors Xi and then averaged. Because this calculation compares only disjoint templates, it will not suffer from the bias introduced by nonindependent templates. This truly unbiased approach has the potentially severe limitation of reducing the number of possible template matches and enlarging the confidence intervals about the SampEn estimate. Because this bias appears to be present only for very small N, the disjoint template approach does not appear necessary in usual practice.

Cross-SampEn shows relative consistency where cross-ApEn does not. As noted above, an essential feature of the measures of order is their relative consistency. That is, if one series is more ordered than another, it should have lower values of ApEn and SampEn for all conditions. We can extend this idea to cross-ApEn and cross-SampEn; if a pair of series is more synchronous than another pair, it should have lower values of cross-ApEn and cross-SampEn statistics for all conditions tested.

We tested the ability of cross-ApEn and cross-SampEn to distinguish between MIX(0.1) and the less-ordered MIX(0.6) processes. The strategy was to compare each with the intermediate MIX(0.3) process. The expected result was that the [MIX(0.3), MIX(0.1)] pair should appear more ordered than the [MIX(0.3), MIX(0.6)] pair, because MIX(0.1) is significantly more ordered than MIX(0.6). That is, cross-ApEn(2, r, 250) [MIX(0.3)∥MIX(0.1)] should be less than cross-ApEn(2, r, 250) [MIX(0.3)∥MIX(0.6)], and cross-SampEn(2, r, 250) [MIX(0.3)∥MIX(0.1)] should be less than cross-SampEn(2, r, 250) [MIX(0.3)∥MIX(0.6)].

We tested this prediction bidirectionally, that is, MIX(0.3) served as the template series for one analysis and as the target series for the next, and over a range of tolerances r by using the bias-max and bias-0 strategies for ensuring that cross-ApEn was defined. The results are shown in Fig. 4. MIX(0.3) was the template series in Fig. 4, C and E, and the target series in Fig. 4, D and F. 


View larger version (28K):
[in this window]
[in a new window]
 
Fig. 4.   Results of cross-ApEn analysis of MIX(P) processes are dep