Estimating the entropy of binary time series: Methodology, some theory and a simulation study

Partly motivated by entropy-estimation problems in neuroscience, we present a detailed and extensive comparison between some of the most popular and effective entropy estimation methods used in practice: The plug-in method, four different estimators based on the Lempel-Ziv (LZ) family of data compression algorithms, an estimator based on the Context-Tree Weighting (CTW) method, and the renewal entropy estimator. **Methodology. Three new entropy estimators are introduced. For two of the four LZ-based estimators, a bootstrap procedure is described for evaluating their standard error, and a practical rule of thumb is heuristically derived for selecting the values of their parameters. ** Theory. We prove that, unlike their earlier versions, the two new LZ-based estimators are consistent for every finite-valued, stationary and ergodic process. An effective method is derived for the accurate approximation of the entropy rate of a finite-state HMM with known distribution. Heuristic calculations are presented and approximate formulas are derived for evaluating the bias and the standard error of each estimator. ** Simulation. All estimators are applied to a wide range of data generated by numerous different processes with varying degrees of dependence and memory. Some conclusions drawn from these experiments include: (i) For all estimators considered, the main source of error is the bias. (ii) The CTW method is repeatedly and consistently seen to provide the most accurate results. (iii) The performance of the LZ-based estimators is often comparable to that of the plug-in method. (iv) The main drawback of the plug-in method is its computational inefficiency.


Introduction
The problem of estimating the entropy of a sequence of discrete observations has received a lot of attention over the past two decades. A, necessarily incomplete, sample of the theory that has been developed can be found in [27] [28].
Information-theoretic methods have been particularly widely used in neuroscience, in a broad effort to analyze and understand the fundamental information-processing tasks performed by the brain. In many of these these studies, the entropy is adopted as the main measure for describing the amount of information transmitted between neurons. There, one of the most basic tasks is to identify appropriate methods for quantifying this information, in other words, to estimate the entropy of spike trains recorded from live animals; see, e.g., [42][51] [43] [35] [26] [18] [30][4] [28] [41].
Motivated, in part, by the application of entropy estimation techniques to such neuronal data, we present a systematic and extensive comparison, both in theory and via simulation, between several of the most commonly used entropy estimation techniques. This work serves, partly, as a more theoretical companion to the experimental work and results presented in [13][12] [14]. There, entropy estimators were applied to the spike trains of 28 neurons recorded simultaneously for a one-hour period from the primary motor and dorsal premotor cortices (MI, PMd) of a monkey. The purpose of those experiments was to examine the effectiveness of the entropy as a statistic, and its utility in revealing some of the underlying structural and statistical characteristics of the spike trains. In contrast, our main aim here is to examine the performance of several of the most effective entropy estimators, and to establish fundamental properties for their applicability, such as rigorous estimates for their convergence rates, bias and variance. In particular, since (discretized) spike trains are typically represented as binary sequences [36], some of our theoretical results and all of our simulation experiments are focused on binary data. Section 2 begins with a description of the entropy estimators we consider. The simplest one is the plug-in or maximum likelihood estimator, which consists of first calculating the empirical frequencies of all words of a fixed length in the data, and then computing the entropy of this empirical distribution. For obvious computational reasons, the plug-in is ineffective for word-lengths beyond 10 or 20, and hence it cannot take into account any potential longer-time dependencies in the data.
A popular approach for overcoming this drawback is to consider entropy estimators based on "universal" data compression algorithms, that is, algorithms which are known to achieve a compression ratio equal to the entropy, for data generated by processes which may possess arbitrarily long memory, and without any prior knowledge about the distri-bution of the underlying process. Since, when trying to estimate the entropy, the actual compression task is irrelevant, many entropy estimators have been developed as modifications of practical compression schemes. Section 2.3 describes two well-known such entropy estimators [22], which are based on the Lempel-Ziv (LZ) family of data compression algorithms [59] [60]. These estimators are known to be consistent (i.e., to converge to the correct value for the entropy) only under certain restrictive conditions on the data. We introduce two new LZ-based estimators, and we prove in Theorem 2.1 that, unlike the estimators in [22], they are consistent under essentially minimal conditions, that is, for data generated by any stationary and ergodic process. Section 2.4 contains an analysis, partly rigorous and partly in terms of heuristic computations, of the rate at which the bias and the variance of each of the four LZ-based estimators converges to zero. A bootstrap procedure is developed for empirically estimating the standard error of two of the four LZ-based estimators, and a practical rule-of-thumb is derived for selecting the values of the parameters of these estimators in practice.
Next, in Section 2.5 we consider an entropy estimator based on the Context-Tree Weighting (CTW) algorithm [53][54] [52]. [In the neuroscience literature, a similar proceedure has been applied in [18] and [26].] The CTW, also originally developed for data compression, can be interpreted as a Bayesian estimation procedure. After a brief description, we explain that it is consistent for data generated by any stationary and ergodic process and show that its bias and variance are, in a sense, as small as can be. 1 Section 3 contains the results of an extensive simulation study, where the various entropy estimators are applied to data generated from numerous different types of processes, with varying degrees of dependence and memory. In Section 3.1, after giving brief descriptions of all these data models, we present (Proposition 3.1) a method for accurately approximating the entropy rate of a Hidden Markov Model (HMM); recall that HMMs are not known to admit closed-form expressions for their entropy rate. Also, again partly motivated by neuroscience applications, we introduce another new entropy estimator, the renewal entropy estimator, which is tailored to binary data generated by renewal processes. Section 3.2 contains a detailed examination of the bias and variance of the four LZbased estimators and the CTW algorithm. There, our earlier theoretical predictions are largely confirmed. Moreover, it is found that two of the four LZ-based estimators are consistently more accurate than the other two.
Finally, Section 3.3 contains a systematic comparison of the performance of all of the above estimators on different types of simulated data. Incidental comparisons between some of these methods on limited data sets have appeared in various places in the literature, and questions have often been raised regarding their relative merits. One of the main goals of the work we report here is to offer clear resolutions for many of these issues. Specifically, in addition to the points mentioned up to now, some of the main conclusion that we draw from the simulation results of Section 3 (these and more are collected in Section 4) can be summarized as follows: • Due its computational inefficiency, the plug-in estimator is the least reliable method, in contrast to the LZ-based estimators and the CTW, which naturally incorporate dependencies in the data at much larger time scales.
• The most effective estimator is the CTW method. Moreover, for the CTW as well as for all other estimators, the main source of error is the bias and not the variance.
• Among the four LZ-based estimators, the two most efficient ones are those with increasing window sizes,Ĥ n of [22] andH n introduced in Section 2.3. Somewhat surprisingly, in several of the simulations we conducted the performance of the LZbased estimators appears to be very similar to that of the plug-in method.

Entropy Estimators and Their Properties
This section contains a detailed description of the entropy estimators that will be applied to simulated data in Section 3. After some basic definitions and notation in Section 2.1, the following four subsections contain the definitions of the estimators together with a discussion of their statistical properties, including conditions for consistency, and estimates of their bias and variance.

Entropy and Entropy Rate
Let X be a random variable or random vector, taking values in an arbitrary finite set A, its alphabet, and with distribution p(x) = Pr{X = x} for x ∈ A. The entropy of X [8] is defined as, where, throughout the paper, log denotes the logarithm to base 2, log 2 . A random process The entropy rate H = H(X), or "per-symbol" entropy, of X is the asymptotic rate at which the entropy of X n 1 changes with n, whenever the limit exists, where H(X 1 , X 2 , . . . , X n ) is the entropy of the jointly distributed random variables X n 1 = (X 1 , X 2 , . . . , X n ). Recall [8] that for a stationary process (i.e., a process such that the distribution of every finite block X n+k n+1 of size k has the same distribution, say p k , independently of its position n), the entropy rate exists and equals, where the conditional entropy H(X n | X n−1 1 ) = H(X n | X n−1 , . . . , X 2 , X 1 ) is defined as, .
As mentioned in the introduction, much of the rest of the paper will be devoted to binary data produced by processes X with alphabet A = {0, 1}.

The Plug-in Estimator
Perhaps the simplest and most straightforward estimator for the entropy rate is the socalled plug-in estimator. Given a data sequence x n 1 of length n, and an arbitrary string, or "word," y w 1 ∈ A w of length w < n, letp w (y w 1 ) denote the empirical probability of the word y w 1 in x n 1 ; that is,p w (y w 1 ) is the frequency with which y w 1 appears in x n 1 . If the data are produced from a stationary and ergodic 2 process, then the law of large numbers guarantees that, for fixed w and large n, the empirical distributionp w will be close to the true distribution p w , and therefore a natural estimator for the entropy rate based on (1) is:Ĥ This is the plug-in estimator with word-length w. Since the empirical distribution is also the maximum likelihood estimate of the true distribution, this is also often referred to as the maximum-likelihood entropy estimator. Suppose the process X is stationary and ergodic. Then, taking w large enough for 1 w H(X w 1 ) to be acceptably close to H in (1), and assuming the number of samples n is much larger than w so that the empirical distribution of order w is close to the true distribution, the plug-in estimatorĤ n,w,plug−in will produce an accurate estimate for the entropy rate. But, among other difficulties, in practice this leads to enormous computational problems because the number of all possible words of length w grows exponentially with w. For example, even for the simple case of binary data with a modest word-length of w = 30, the number of possible strings y w 1 is 2 30 , which in practice means that the estimator would either require astronomical amounts of data to estimatep w accurately, or it would be severely undersampled. See also [31] and the references therein for a discussion of the undersampling problem. Another drawback of the plug-in estimator is that it is hard to quantify its bias and to correct for it. For any fixed word-length w, it is easy to see that the bias E Ĥ n,w,plug−in − 1 w H(X w 1 ) is always negative [1] [30], whereas the difference between the wth-order persymbol entropy and the entropy rate, 1

The Lempel-Ziv Estimators
An intuitively appealing and popular way of estimating the entropy of discrete data with possibly long memory, is based on the use of so-called universal data compression algorithms. These are algorithms that are known to be able to optimally compress data from an arbitrary process (assuming some broad conditions are satisfied), where optimality means that the compression ratio they achieve is asymptotically equal to the entropy rate of the underlying process -although the statistics of this process are not assumed to be known a priori. Perhaps the most commonly used methods in this context are based on a family of compression schemes known as Lempel-Ziv (LZ) algorithms; see, e.g., [59][60] [56].
Since the entropy estimation task is simpler than that of actually compressing the data, several modified versions of the original compression algorithms have been proposed and used extensively in practice. All these methods are based on the calculation of the lengths of certain repeating patterns in the data. Specifically, given a data realization x = (. . . , x −1 , x 0 , x 1 , x 2 , . . .), for every position i in x and any "window length" n ≥ 1, consider the length ℓ of the longest segment x i+ℓ−1 i in the data starting at i which also appears in the window x i−1 i−n of length n preceding position i. Formally, define L n i = L n i (x) = L n i (x i+n−1 i−n ) as 1+ [that longest match-length]: Suppose that the process X is stationary and ergodic, and consider the random match- . In [56] [29] it was shown that, for any fixed position i, the match-lengths grow logarithmically with the window size n, and in fact, where H is the entropy rate of the process. This result suggests that the quantity (log n)/L n i can be used as an entropy estimator, and, clearly, in order to make more efficient use of the data and reduce the variance, it would be more reasonable to look at the average value of various match-lengths L n i taken at different positions i; see the discussion in [22]. To that effect, the following two estimators are considered in [22]. Given a data realization x = x ∞ −∞ , a window length n ≥ 1, and a number of matches k ≥ 1, the sliding-window LZ estimatorĤ n,k =Ĥ n,k (x) =Ĥ n,k (x n+k−1 −n+1 ) is defined by, Similarly, the increasing-window LZ estimatorĤ n =Ĥ n (x) =Ĥ n (x 2n−1 0 ) is defined by, The difference between the two estimators in (3) and (4) is thatĤ n,k uses a fixed window length, whileĤ n uses the entire history as its window, so that the window length increases as the matching position moves forward.
In [22] it is shown that, under appropriate conditions, both estimatorsĤ n,k andĤ n are consistent, in that they converge to the entropy rate of the underlying process with probability 1, as n, k → ∞. Specifically, it is assumed that the process is stationary and ergodic, that it takes on only finitely many values, and that it satisfies the Doeblin Condition (DC). This condition says that there is a finite number of steps, say r, in the process, such that, after r time steps, no matter what has occurred before, anything can happen with positive probability: Doeblin Condition (DC). There exists an integer r ≥ 1 and a real number β > 0 such that, Pr(X r = a | X 0 −∞ ) > β, for all a ∈ A and with probability one in the conditioning, i.e., for almost all semi-infinite realizations of the past X 0 −∞ = ( X 0 , X −1 , . . .).
Condition (DC) has the advantage that it is not quantitative -the values of r and β can be arbitrary -and, therefore, for specific applications it is fairly easy to see whether it is satisfied or not. But it is restrictive, and, as it turns out, it can be avoided altogether if we consider a modified version of the above two estimators.
To that end, we define two new estimatorsH n,k andH n as follows. Given x = x ∞ −∞ , n and k as above, define the new sliding-window estimatorH n,k =H n,k (x) =H n,k (x n+k−1 −n+1 ), and the new increasing-window estimatorH n =H n (x) =H n (x 2n−1 0 ) as, Below some basic properties of these four estimators are established, and conditions are given for their asymptotic consistency. Parts (i) and (iii) of Theorem 2.1 are new; most of part (ii) is contained in [22]. (i) When applied to an arbitrary data string, the estimators defined in (3)-(6) always satisfy,Ĥ n,k ≤H n,k andĤ n ≤H n , for any n, k.
(ii) The estimatorsĤ n,k andĤ n are consistent when applied to data generated by a finite-valued, stationary, ergodic process that satisfies Doeblin's condition (DC). With probability one we have:Ĥ n,k → H,Ĥ n → H, as k, n → ∞.
(iii) The estimatorsH n,k andH n are consistent when applied to data generated by an arbitrary finite-valued, stationary, ergodic process, even if (DC) does not hold. With probability one we have:H n,k → H,H n → H, as k, n → ∞.
Note that parts (ii) and (iii) do not specify the manner in which n and k go to infinity. The results are actually valid in the following cases: 1. If the two limits as n and k tend to infinity are taken separately, i.e., first k → ∞ and then n → ∞, or vice versa; 2. If k → ∞ and n = n k varies with k in such a way that n k → ∞ as k → ∞; 3. If n → ∞ and k = k n varies with n in such a way that it increases to infinity as n → ∞; 4. If k and n both vary arbitrarily in such a way that k stays bounded and n → ∞.
Proof. Part (i). An application of Jensen's inequality to the convex function x → 1/x, with respect to the uniform distribution (1/k, . . . , 1/k) on the set {1, 2, . . . , k}, yields, as required. The proof of the second assertion is similar.
Part (ii). The results here are, for the most part, proved in [22], where it is established thatĤ n → H andĤ n,n → H as n → ∞, with probability one. So it remains to show that H n,k → H as n, k → ∞ in each of the four cases stated above. For case 1 observe that, with probability 1, where (a) follows from (2). To reverse the limits, we define, for each fixed n, a new process {Y i } is itself stationary and ergodic. Recalling also from [22] that the convergence in (2) takes place not only with probability one but also in L 1 , we may apply the ergodic theorem to obtain that, with probability 1, where (b) follows by the ergodic theorem and (c) from the L 1 version of (2). The proof of case 2 is identical to the case k = n considered in [22]. In case 3, since the sequence {k n } is increasing, the limit of H kn,n reduces to a subsequence of the corresponding limit in case 2 upon considering the inverse sequence {n k }.
Finally for case 4, recall from (7) that, lim n H k,n = H with prob. 1, for any fixed k. Therefore the same will hold with a varying k, as long as it varies among finitely many values.
Part (iii). The proofs of the consistency results forH n andH n,k can be carried out along the same lines as the proofs of the corresponding results in [22], together with their extensions as in Part (ii) above. The only difference is in the main technical step, namely, the verification of a uniform integrability condition. In the present case, what is needed is to show that, This is done in the following lemma. Proof. Given a data realization x = (. . . , x −2 , x −1 , x 0 , x 1 , x 2 , . . .) and an m ≥ 1, the recurrence time R m is defined as the first time the substring x m 1 appears again in the past. More precisely, R m is the number of steps to the left of x m 1 we have to look in order to find a copy of x m 1 : For any such realization x and any n ≥ 1, if we take m = L n 1 , then by the definitions of R m and L n 1 it follows that, R m > n, which implies, log R m m > log n L n 1 , and thus it is always the case that, Therefore, to establish (8) it suffices to prove: To that end, we expand this expectation as, where K is an arbitrary integer to be chosen later. Applying Markov's inequality, To calculate the expectation of R m , suppose that the process X takes on α = |A| possible values, so that there are α m possible strings x m 1 of length m. Now recall Kac's theorem [56] Combining (10) and (11) yields, where we choose K > log α. This establishes (9) and completes the proof. 2

Bias and Variance of the LZ-based Estimators
In practice, when applying the sliding-window LZ estimatorsĤ n,k orH n,k on finite data strings, the values of the parameters k and n need to be chosen, so that k + n is approximately equal to the total data length. This presents the following dilemma: Using a long window size n, the estimators are more likely to capture the longer-term trends in the data, but, as shown in [33][45], the match-lengths L n i starting at different positions i have large fluctuations. So a large window size n and a small number of matching positions k will give estimates with high variance. On the other hand, if a relatively small value for n is chosen and the estimate is an average over a large number of positions k, then the variance will be reduced at the cost of increasing the bias, since the expected value of L n i / log n is known to converge to 1/H very slowly [55]. Therefore n and k need to be chosen in a way such that the above bias/variance tradeoff is balanced. From the earlier theoretical results of [33][45] [46][55] [57] it follows that, under appropriate conditions, the bias is approximately of the order O(1/ log n), whereas from the central limit theorem it is easily seen that the variance is approximately of order O(1/k). This indicates that the relative values of n and k should probably be chosen to satisfy k ≈ O((log n) 2 ).
Although this is a useful general guideline, we also consider the problem of empirically evaluating the relative estimation error on particular data sets. Next we outline a bootstrap procedure, which gives empirical estimates of the varianceĤ n,k andH n,k ; an analogous method was used for the estimatorĤ n,k in [44], in the context of estimating the entropy of whale songs.
Let L denote the sequence of match-lengths L = (L n 1 , , L n 2 , . . . , L n k ) computed directly from the data, as in the definitions ofĤ n,k andH n,k . Roughly speaking, the proposed procedure is carried out in three steps: First, we sample with replacement from L in order to obtain many pseudo-time series with the same length as L; then we compute new entropy estimates from each of the new sequences usingĤ n,k orH n,k ; and finally we estimate the variance of the initial entropy estimates as the sample variance of the new estimates. The most important step is the sampling, since the elements of sequence (L n 1 , . . . , L n k ) are not independent. In order to maintain the right form of dependence, we adopt a version of the stationary bootstrap procedure of [34]. The basic idea is, instead of sampling individual L n i 's from L, to sample whole blocks with random lengths. The choice of the distribution of their lengths is made in such a way as to guarantee that they are typically long enough to maintain sufficient dependence as in the original sequence. The results in [34] provide conditions which justify the application of this procedure.
The details of the three steps above are as follows: First, a random position j ∈ {1, 2, . . . , k} is selected uniformly at random, and a random length T is chosen with geometric distribution with mean 1/p (the choice of p is discussed below). Then the block of match-lengths (L n j , L n j+1 , . . . , L j+T −1 ) is copied from L, and the same process is repeated until the concatenation L * of the sampled blocks has length k. This gives the first bootstrap sample. Then the whole process is repeated to generate a total of B such blocks L * 1 , L * 2 , . . . , L * B , each of length k. From these we calculate new entropy estimatesĤ * (m) orH * (m), for m = 1, 2, . . . , B, according to the definition of the entropy estimator being used, as in (3) or (5), respectively; the choice of the number B of blocks is discussed below. The bootstrap estimate of the variance ofĤ is, m=1Ĥ * (m); similarly forH. The choice of the parameter p depends on the length of the memory of the matchlength sequence L = (L n 1 , L n 2 , . . . , L n k ); the longer the memory, the larger the blocks need to be, therefore, the smaller the parameter p. In practice, p is chosen by studying the autocorrelogram of L, which is typically decreasing with the lag: We choose an appropriate cutoff threshold, take the corresponding lag to be the average block size, and choose p as the reciprocal of that lag. Finally, the number of blocks B is customarily chosen large enough so that the histogram of the bootstrap samplesĤ * (1),Ĥ * (2), . . . ,Ĥ * (B) "looks" approximately Gaussian. Typical values used in applications are between B = 500 and B = 1000. In all our experiments in [12] [14] and in the results presented in the following section we set B = 1000, which, as discussed below, appears to have been sufficiently large for the central limit theorem to apply to within a close approximation.

Context-Tree Weighting
One of the fundamental ways in which the entropy rate arises as a natural quantity, is in the Shannon-McMillan-Breiman theorem [8]; it states that, for any stationary and ergodic process X = {. . . , X −1 , X 0 , X 1 , X 2 , . . .} with entropy rate H, where p n (X n 1 ) denotes the (random) probability of the random string X n 1 . This suggests that, one way to estimate H from a long realization x n 1 of X, is to first somehow estimate its probability and then use the estimated probabilityp n (x n 1 ) to obtain an estimate for the entropy rate via,Ĥ The Context-Tree Weighting (CTW) algorithm [53][54] [52] is a method, originally developed in the context of data compression, which can be interpreted as an implementation of hierarchical Bayesian procedure for estimating the probability of a string generated by a binary "tree process." 3 The details of the precise way in which the CTW operates can be found in [53][54] [52]; here we simply give a brief overview of what (and not how) the CTW actually computes. In order to do that, we first need to describe tree processes.
A binary tree process of depth D is a binary process X with a distribution defined in terms of a suffix set S, consisting of binary strings of length no longer than D, and a parameter vector θ = (θ s ; s ∈ S), where each θ s ∈ [0, 1]. The suffix set S is assumed to be complete and proper, which means that any semi-infinite binary string x 0 −∞ has exactly one suffix s is S, i.e., there exists exactly one s ∈ S such that x 0 −∞ can be written as x −k −∞ s, for some integer k. 4 We write s = s(x 0 −∞ ) ∈ S for this unique suffix. Then the distribution of X is specified by defining the conditional probabilities, . It is clear that the process just defined could be thought of simply as a D-th order Markov chain, but this would ignore the important information contained in S: If a suffix string s ∈ S has length ℓ < D, then, conditional on any past sequence x n −∞ which ends in s, the distribution of X n+1 only depends on the most recent ℓ symbols. Therefore, the suffix set offers an economical way for describing the transition probabilities of X, especially for chains that can be represented with a relatively small suffix set.
Suppose that a certain string x n 1 has been generated by a tree process of depth no greater than D, but with unknown suffix set S * and parameter vector θ * = (θ * s ; s ∈ S). Following classical Bayesian methodology, we assign a prior probability π(S) on each (complete and proper) suffix set S of depth D or less, and, given S, we assign a prior probability π(θ | S) on each parameter vector θ = (θ s ). A Bayesian approximation to the true probability of x n 1 (under S * and θ * ) is the mixture probability, where P S,θ (x n 1 ) is the probability of x n 1 under the distribution of a tree process with suffix set S and parameter vector θ. The expression in (14) is, in practice, impossible to compute directly, since the number of suffix sets of depth ≤ D (i.e., the number of terms in the sum) is of order 2 D . This is obviously prohibitively large for any D beyond 20 or 30.
The CTW algorithm is an efficient procedure for computing the mixture probability in (14), for a specific choice of the prior distributions π(S), π(θ | S): The prior on S is, where |S| is the number of elements of S and N (S) is the number of string is S with length strictly smaller than D. Given a suffix set S, the prior on θ is the product ( 1 2 , 1 2 )-Dirichlet distribution, i.e., under π(θ|S) the individual θ s are independent, with each θ s ∼ Dirichlet( 1 2 , 1 2 ). The main practical advantage of the CTW algorithm is that it can actually compute the probability in (14) exactly. In fact, this computation can be performed sequentially, in linear time in the length of the string n, and using an amount of memory which also grows linearly with n. This, in particular, makes it possible to consider much longer memory lengths D than would be possible with the plug-in method.
The CTW Entropy Estimator. Thus motivated, given a binary string x n 1 , we define the CTW entropy estimatorĤ n,D, ctw as, whereP D, mix (x n 1 ) is the mixture probability in (14) computed by the CTW algorithm. [Corresponding proceedures are similarly described in [18] [26].] The justification for this definition comes from the discussion leading to equation (13) above. Clearly, if the true probability of x n 1 is P * (x n 1 ), the estimator performs well when, In many cases this approximation can be rigorously justified, and in certain cases it can actually be accurately quantified.
Assume that x n 1 is generated by an unknown tree process (of depth no greater than D) with suffix set S * . The main theoretical result of [53] states that, for any string x n 1 , of any finite length n, generated by any such process, the difference between the two terms in (16) can be uniformly bounded above; from this it easily follows that, This nonasymptotic, quantitative bound, easily leads to various properties ofĤ n,D, ctw : First, (17) combined with the Shannon-McMillan-Breiman theorem (12) and the pointwise converse source coding theorem [2][19], readily implies thatĤ N,D, ctw is consistent, that is, it converges to the true entropy rate of the underlying process, with probability one, as n → ∞. Also, Shannon's source coding theorem [8] implies that the expected value ofĤ n,D, ctw cannot be smaller than the true entropy rate H; therefore, taking expectations in (17) gives, In view of Rissanen's [37] well-known universal lower bound, (18) shows that the bias of the CTW is essentially as small as can be. Finally, for the variance, if we subtract H from both sides of (17), multiply by √ n, and apply the central-limit refinement to the Shannon-McMillan-Breiman theorem [58][ 16], we obtain that the standard deviation of the estimatesĤ n,D, ctw is ≈ σ X / √ n, where σ 2 X is the minimal coding variance of X [21]. This is also optimal, in view of the second-order coding theorem of [21].
Therefore, for data x n 1 generated by tree processes, the bias of the CTW estimator is of order O(log n/n), and its variance is O(1/n). Compared to the earlier LZ-based estimators, these bounds suggest much faster convergence, and are in fact optimal. In particular, the O(log n/n) bound on the bias indicates that, unlike the LZ-based estimators, the CTW can give useful results even on small data sets.
Extensions. An important issue for the performance of the CTW entropy estimator, especially when used on data with potentially long-range dependence, is the choice of the depth D: While larger values of D give the estimator a chance to capture longer-term trends, we then pay a price in the algorithm's complexity and in the estimation bias. This issue will not be discussed further here; a more detailed discussion of this point along with experimental results can be found in [12] [14].
The CTW algorithm has also been extended beyond finite-memory processes [52]. The basic method is modified to produce an estimated probabilityP ∞ (x n 1 ), without assuming a predetermined maximal suffix depth D. The sequential nature of the computation remains exactly the same, leading to a corresponding entropy estimator defined analogously to the one in (15), asĤ n,∞, ctw = − 1 n logP ∞ (x n 1 ). Again it is easy to show thatĤ n,∞, ctw is consistent with probability one, this time for every stationary and ergodic (binary) process. The price of this generalization is that the earlier estimates for the bias and variance no longer apply, although they do remain valid is the data actually come from a finite-memory process. In numerous simulation experiments we found that there is no significant advantage in usingĤ n,D, ctw with a finite depth D overĤ n,∞, ctw , except for the somewhat shorter computation time. For that reason, in all the experimental results in Sections 3.2 and 3.3 below, we only report estimates obtained byĤ n,∞, ctw .
Finally, perhaps the most striking feature of the CTW algorithm is that it can be modified to compute the "best" suffix set S that can be fitted to a given data string, where, following standard statistical (Bayesian) practice, "best" here means the one which is most likely under the posterior distribution. To be precise, recall the prior distributions π(S) and π(θ | S) on suffix sets S and on parameter vectors θ, respectively. Using Bayes' rule, the posterior distribution on suffix sets S can be expressed, ; the suffix setŜ =Ŝ(x n 1 ) which maximizes this probability is called the Maximum A posteriori Probability, or MAP, suffix set. Although the exact computation ofŜ is, typically, prohibitively hard to carry out directly, the Context-Tree Maximizing (CTM) algorithm proposed in [50] is an efficient procedure (with complexity and memory requirements essentially identical to the CTW algorithm) for computingŜ. The CTM algorithm will not be used or discussed further in this work; see the discussion in [12] [14], where it plays an important part in the analysis of neuronal data.

Results on Simulated Data
This section contains the results of an extensive simulation study, comparing various aspects of the behavior of the different entropy estimators presented earlier (the plug-in estimator, the four LZ-based estimators, and the CTW estimator), applied to different kinds of simulated binary data. Motivated, in part, by applications in neuroscience, we also introduce a new method, the renewal entropy estimator. Section 3.1 contains descriptions of the statistical models used to generate the data, along with exact formulas or close approximations for their entropy rates. Sections 3.2 and 3.3 contain the actual simulation results.

I.I.D. (or "homogeneous Poisson") Data
The simplest model of a binary random process X is as a sequence of independent and identically distributed (i.i.d.) random variables {X n }, where the the X i are independent and all have the same distribution. This process has no memory, and in the neuroscience literature it is often referred to as a "homogeneous Poisson process." The distribution of each X i is described by a parameter p ∈ [0, 1], so that Pr{X i = 1} = p and Pr{X i = 0} = 1 − p. If {X n } were to represent a spike train, then p = E(X i ) would be its average firing rate.
The entropy rate of this process is simply,

Markov Chains
An ℓ-th order (homogeneous) Markov chain with (finite) alphabet A is a process X = {X n } with the property that, for all x n 1 ∈ A n , where P = P x ℓ 1 ,x ℓ+1 is the transition matrix of the chain. This formalizes the idea that the memory of the process has length ℓ: The probability of each new symbol depends only on the most recent ℓ symbols, and it is conditionally independent of the more distant past. The distribution of this process is described by the initial distribution of its first ℓ symbols, (π(x ℓ 1 )), and the transition matrix P . The entropy rate of an ergodic (i.e., irreducible and aperiodic) ℓ-th order Markov chain X is given by, where π * is the unique stationary distribution of the chain.

Hidden Markov Models
Next we consider a class of binary processes X = {X n }, called hidden Markov models (HMMs) or hidden Markov processes, which typically have infinite memory. For the purposes of this discussion, a binary HMM can be defined as follows. Suppose Y = {Y n } is a first-order, ergodic Markov chain, which is stationary, i.e., its initial distribution π is the same as its stationary distribution π * . Let the alphabet of Y be an arbitrary finite set A, write P = (P yy ′ ) for its transition matrix, and let Q = (Q yx ; y ∈ A, x ∈ {0, 1}) be a different transition matrix, from A to {0, 1}. Then, for each n, given {Y i } and the previous values x n−1 1 of X n−1

1
, the distribution of the random variable X n is, independently of the remaining {Y i } and of X n−1

1
. The resulting process X = {X n } is a binary HMM.
The consideration of HMMs here is partly motivated by the desire to simulate spike trains with slowly varying rates, as in the case of real neuronal firing. To illustrate, consider the following (somewhat oversimplified) description of a model that will be used in the simulation examples below. Let the Markov chain Y = {Y n } represent the process which modulates the firing rate of the binary "spike train" X, so that Y takes a finite number of values, A = {r 1 , r 2 , . . . , r α }, with each r i ∈ (0, 1). These values correspond to α different firing regimes, so that, e.g., Y n = r 1 means that the average firing rate at that instant is r 1 spikes-per-time-unit. To ensure that the firing rate varies slowly, define, for every y ∈ A, the transition probability that Y n remains in the same state to be Pr{Y n = r | Y n−1 = r} = 1 − ǫ, for some small ǫ > 0. Then, conditional on {Y n }, the distribution of each X n is given by Pr{X n = 1 | Y n = r} = r = 1 − Pr{X n = 0 | Y n = y}.
In general, an HMM defined as above is stationary, ergodic and typically has infinite memory -it is not a ℓ-th order Markov chain for any ℓ; see [10] and the references therein for details. Moreover, there is no closed-form expression for the entropy rate of a general HMM, but, as outlined below, it is fairly easy to obtain an accurate approximation when the distribution of the HMM is known a priori, via the Shannon-McMillan-Breiman theorem (12). That is, the value of the entropy rate H = H(X) can be estimated accurately as long as it is possible to get a close approximation for the probability p n (X n 1 ) of a long random sample X n 1 . This calculation is, in principle, hard to perform, since it requires the computation of an average over all possible state sequences y n 1 , and their number grows exponentially with n. As it turns out, adapting an idea similar to the usual dynamic programming algorithm used for HMM state estimation, the required probability p n (X n 1 ) can actually be computed very efficiently; similar techniques appear in various places in the literature, e.g., [10] [17]. Here we adopt the following method, developed independently in [12].
First, generate and fix a long random realization x n 1 of the HMM X, and define the matrices M (k) by, The following proposition says that the probability of an arbitrary x n 1 can be obtained in n matrix multiplications, involving matrices of dimension |A| × |A|. For moderate alphabet sizes, this can be easily carried out even for large n, e.g., on the order of 10 6 . Moreover, as the HMM process X inherits the strong mixing properties of the underlying Markov chain Y , Ibragimov's central-limit refinement [16] to the Shannon-McMillan-Breiman theorem suggests that the variance of the estimates obtained should decay at the rapid rate of 1/n. Therefore, in practice, it should be possible to efficiently obtain a reliable, stable approximation for the entropy rate, a prediction we have repeatedly verified through simulation. 5 Proposition 3.1. Under the above assumptions, the probability of an arbitrary x n 1 ∈ {0, 1} n produced by the HMM X can be expressed as, where 1 is the column vector of |A| 1s.
Proof: Let y n 1 ∈ A n denote any specific realization of the hidden Markov process Y n 1 . We have,

Renewal Processes
A common alternative mathematical description for the distribution of a binary string is via the distribution of the time intervals between successive 1s, or, in the case of a binary spike train, the interspike intervals (ISIs) as they are often called in the neuroscience literature [42] [4]. Specifically, let {t i } denote the sequence of times t when X t = 1, and let {Y i = t i+1 − t i } be the sequence of "interarrival times" or ISIs of X = {X n }. Instead of defining the joint distributions of blocks X n 1 of random variables from X, its distribution can be specified by that of the process Y = {Y i }. For example, if the Y i are i.i.d. random variables with geometric distribution with parameter p ∈ [0, 1], then X is itself a (binary) i.i.d. process with parameter p.
More generally, a renewal process X is defined in terms of an arbitrary i.i.d. ISI process Y , with each Y i having a common discrete distribution P = (p j ; j = 1, 2, . . .). 6 Recall [32] that the entropy rates H(X) of X and H(Y ) of Y are related by, where λ = E(X 1 ) = 1/E(Y 1 ). This simple relation motivates the consideration of a different estimator for the entropy rate of a renewal process. Given binary data x n 1 of length n, calculate the corresponding sequence of ISIs y m 1 , and define, is the empirical distribution of the ISIs y m 1 , andλ is the empirical rate of X, i.e., the proportion of 1s in x n 1 . We call this the renewal entropy rate estimator. When the data are indeed generated by a renewal process, the law of large numbers guarantees thatĤ n,renewal will converge to H(X), with probability one, as n → ∞. Moreover, under rather weak regularity conditions on the distribution of the ISIs, the central limit theorem implies that the variance of these estimates decays like O(1/n). But the results of the renewal entropy estimator remain meaningful even if X is not a renewal process. For example, if the corresponding ISI process Y is stationary and ergodic but not i.i.d., then the renewal entropy estimator will converge to the value λH(Y 1 ), whereas the true entropy rate in this case is the (smaller) quantity, Therefore, the renewal entropy estimator can be employed to test for the presence or absence of renewal structure in particular data sets, by comparing the value ofĤ n,renewal with that of other estimators; see [14][12] for a detailed such study.

Bias and Variance of the CTW and the LZ-based Estimators
In order to compare the empirical bias and variance of the four LZ-based estimators and the CTW estimator with the corresponding theoretical predictions of Sections 2.4 and 2.5, the five estimators are applied to simulated (binary) data. [In this as well as in the following section, we only present results for the infinite-suffix-depth CTW estimatorĤ n,∞, ctw ; recall the discussion at the very end Section 2.5.] The simulated data were generated from four different processes: An i.i.d. (or "homogeneous Poisson") process, and three different Markov chains with orders ℓ = 1, 2, and 10. For each parameter setting of each model, 50 independent realizations were generated and the bias of each method was estimated by subtracting the true entropy rate from the empirical mean of the individual estimates. Similarly, the standard error was estimated by the empirical standard deviation of these results.
Specifically, forĤ n,k andH n,k , first a window size n and a number of matches k were selected, and then 50 independent realizations of length N = 2n + k were generated. The bias and standard error are plotted in Figure 1 against O(1/ log n) and O(1/ √ k), respectively. The approximately linear curves confirm the theoretical predictions that the bias and variance decay to zero at rates O(1/ log n) and O(1/k), respectively. Note that, althoughH n,k is always larger thanĤ n,k , there is no systematic trend regarding which one gives more accurate estimates. Also plotted on Figure 1 is the bias of infinite-suffix-depth CTW estimator,Ĥ N,∞,ctw , applied to data with total length N = 2n + k.  Figure 1: Results obtained byĤ n,k (shown as red lines with circles), byH n,k (blue lines with triangles), and by the CTW estimatorĤ N,∞, ctw (green lines with squares), applied to data from four different processes. ForĤ n,k andH n,k , the first row shows the bias plotted against 1/ log n, for window lengths n = 10 3 , 10 4 , 10 5 , 10 6 , and with a fixed number of matches k = 10 6 ; the second row shows the standard error plotted against 1/ √ k, where k = 10 4 , 10 5 , 10 6 , with fixed n = 10 6 . ForĤ N,∞,ctw , the first row shows the bias when applied to data of the same total length N = n + 2k; for its standard error see Figure 2. The results of the CTW are generally much more accurate than those of the LZ estimators, except in the case of the 10th order Markov chain with small data size, where the LZ-based methods seem to get better results.
The values of the bias and standard error ofĤ n,k andH n,k in Figure 1 suggest that, in order to minimize the total mean squared error (MSE) of the LZ-based estimates, the number of matches k should be chosen to be small relative to the window length n, since the variance clearly appears to decay much faster than the square of the bias. This is further confirmed by the results shown in Table 1, which shows corresponding estimates for an i.i.d. process with p = 0.25 and data length N = 10 6 . The values of n and k satisfy n + k = N − 2 log N . The ratio n/k ranges from 1 to 10,000.
Next, the validity of the bootstrap procedure for the standard error of the LZ estimatorŝ H n,k andH n,k is examined; recall the description in Section 2.4. For data generated from the same four processes, the bootstrap estimate of the standard error is compared to that  Table 2: Comparison between the bootstrap standard error and the empirical estimate of the standard deviation for the four different processes. The window size n is fixed at 10 3 .
Corresponding experiments were performed forĤ n andH n , applied to data from the same four types of processes. Although for these estimators there is little in the way of rigorous theory that can be used as a guideline to compute their mean and variance, simple heuristic calculations strongly suggest that they should decay like O(1/ log n) and O(1/n), respectively. Figure 2 shows the bias and standard error computed empirically as forĤ n,k andH n,k , and plotted against 1/ log n and 1/ √ n, respectively. Once again, the fact that the resulting empirical curves are approximately linear agrees with these predictions. Observe that the bias ofH n can be either positive or negative; the same behavior was observed forĤ n in numerous simulation experiments. The fact that the standard error converges much faster than the bias strongly suggests that, for all four LZ-based estimators, it is the bias that dominates the estimation error,  Table 3: Comparison between the bootstrap standard error and the empirical estimate of the standard deviation for the four different processes. The window size n is fixed at 10 3 .
even in these simple cases of processes with fairly short memory. Figure 2 also shows the bias and standard error of the CTW estimatorĤ n,∞, ctw . Its bias appears to generally converge significantly faster than that of the LZ-based methods, while its standard error is very close to that ofĤ n andH n , apparently converging to zero at a very similar rate. Nevertheless, the results show that the main source of error of the CTW is again from the bias.
Finally we note that, in the results shown, as well as in numerous other simulation runs, we found that the two increasing-window LZ estimatorsĤ n andH n generally performed better, and always at least as well, as the corresponding sliding-window estimatorsĤ n,k andH n,k . Specifically, the simulation results show that the optimal choice of parameters forĤ n,k andH n,k leads to estimates whose bias is very similar to that ofĤ n andH n , whereas the increasing-window estimators have significantly smaller variance.

Comparison of the Different Entropy Estimators
This section contains a systematic comparison of the performance of all the estimators introduced so far: The plug-in method, the LZ-based estimators, the CTW algorithm, and the renewal entropy estimator. All the estimators are applied to different types of simulated binary data, generated from processes with different degrees of dependence and memory. The main figure of merit adopted here is the ratio √ MSE H , between the squareroot of the mean square error (MSE), and the true entropy rate. The corresponding ratios of the bias to the entropy and of the standard error to the entropy are also considered.
Since, as noted earlier, the two increasing-window LZ estimatorsĤ n ,H n generally perform better or at least as well as their sliding-window counterparts,Ĥ n,k ,H n,k , from now on we restrict attention toĤ n andH n . Table 4 shows the results obtained by the plug-in for word-lengths w = 15 and w = 20, the LZ-estimatorsĤ n andH n , and the CTW estimatorĤ n,∞, ctw , on data generated from the same four finite-memory processes as above: An i.i.d. process and three Markov chains of order ℓ = 1, 2 and 10. Again we observe that the main source of error is from the bias for  Figure 2: Results obtained byĤ n (shown as red lines with circles), byH n (blue lines with triangles), and by the CTW estimator (green lines with squares), applied to data from four different processes. ForĤ n andH n , the first row shows the bias plotted against 1/ log n, for window lengths n = 5000, 10 4 , 10 5 , 10 6 , and the second row shows the standard error plotted against 1/ √ n. Similarly, forĤ n,∞,ctw , the two rows show the bias and standard error when the estimator is applied to data of the same total length. The true values of the entropy rates of the four processes are the same as before.
all five estimators. The CTW appears to have the smallest bias, while the standard error varies little among the different estimators. More specifically, the results demonstrate that, for short memory processes, the plug-in estimates are often better than those obtained by the LZ estimators, whereas for the 10-th order Markov chain the plug-in with word-length w = 20 is much worse thanĤ n andH n because of the undersampling problem mentioned earlier. The CTW estimator performs uniformly better than all the other estimators, for both short and relatively long memory processes. Its fast convergence rate outperforms the LZ-based estimators, and its ability to look much further into the past makes it much more accurate than the plug-in. Table 5 shows corresponding results for data generated by two different binary HMMs; recall the relevant description from Section 3.1.3. The first one has three (hidden) states, A = {r 1 , r 2 , r 3 }, where r 1 = 0.005, r 2 = 0.02, r 3 = 0.05. The transition matrix of Y has Pr{Y n+1 = r | Y n = r} = 1 − ǫ, and Pr{Y n+1 = r ′ | Y n = r} = ǫ/2, for any r and r ′ = r, where ǫ = 0.001. Given the sequence Y = {Y n } of hidden states, the observations X = {X n } are conditionally independent samples, where each X n is a binary random variable with Pr{X n = 1|Y n = r} = 1 − Pr{X n = 0|Y n = r} = r. In the second example, the hidden process Y = {Y n } takes values in the set A of 50 real numbers evenly spaced between 0.001 and 0.1. The evolution of the process Y = {Y n } is that of a nearest-  Table 4: Comparison of the ratios between the bias, the standard error and the square-root of the MSE over the true entropy rate. Five estimators are used on four different types of data. All estimators are applied to data of the same length n = 10 6 . Results are shown as percentages, that is, all ratios are multiplied by 100.
neighbor random walk; Y n stays constant with probability (1 − ǫ) and it moves to either one of its two neighbors with probability ǫ/2, with ǫ = 0.02. The conditional distribution of X given Y is the same as before. In both cases, the "true" entropy of the underlying HMM was computed using the approximation procedure described earlier. In the first example, 10 independent realizations of the HMM were used according to the formula of Proposition 3.1, and the average of the resulting estimates was taken to be the true entropy rate in the calculations of the bias and standard error in Table 5. The same procedure was applied to three independent realizations of the second example.
The results on HMM data are quite similar to those obtained on processes with short memory shown in Table 4. The LZ-based estimators have heavy biases, whereas both the plug-in and the CTW method give very accurate estimates. This is probably due to the fact that, although HMMs in general have infinite memory, the memory of an HMM with a finite number of states decays exponentially, therefore only the short-range statistical dependence in the data is significant.
To simulate binary data strings with potentially longer memory (and with characteristics that are generally closer to real neuronal spike trains), we turn to renewal processes X = {X n } with ISIs Y = {Y i } that are distributed according to a mixture of discretized Gamma distributions; recall the description of Section 3.  Table 5: Comparison of the ratios between the bias, the standard error and the squareroot of the MSE over the true entropy rate. Five estimators are used on data generated from two different HMMs. All estimators are applied to data of the same length n = 10 6 . Results are shown as percentages, that is, all ratios are multiplied by 100.
i.i.d. with distribution P = (p j ) given by, where µ ∈ (0, 1) is the mixing proportion and each f i is a (discretized) Gamma density with parameters (α i , β i ), respectively. If the ISI distribution P was simply Gamma(α, β), then the rate E(X 1 ) of the binary process X would be the reciprocal of E(Y 1 ) = αβ; therefore, taking a mixture of two Gamma densities, one with αβ small and one with αβ large, gives an approximate model of a "bursty" neuron, that is, a neuron which sometimes fires several times in rapid succession, and sometimes does not fire for a long period. Figure 3 shows the result of such a simulated process X. The empirical ISI density resembles a real neuron's ISI distribution, and the peak near the zero lag in the autocorrelogram of {X n } indicates the bursty behavior. Table 6 shows the results of the estimates obtained by the renewal entropy estimator, the plug-in with word-length w = 20, the two increasing-window LZ-based estimators, and the CTW method, applied to data generated by a binary renewal process. The ISI distribution P is a mixture of two (discrete) Gamma densities f 1 and f 2 , where f 1 has fixed parameters (α 1 , β 1 ) = (2, 10) that represent the bursting regime, and f 2 takes three different sets of (much larger) parameters (α 2 , β 2 ), representing low frequency firing.
The CTW and renewal estimator consistently outperform the other methods, and, on the other extreme, the LZ-based estimators have high biases and give relatively poor results. The plug-in estimator only considers what happens in 20-bit-long windows, and therefore it misses all the data features beyond this range. Especially in the case of large parameters (α 2 , β 2 ) where the ISI distribution has heavier tails, the structure and the dependence in the data becomes significantly longer in range. Also, in that regime, the resulting number of ISI data points y i is much smaller, and therefore, the empirical ISI distribution is less accurate; as a result, the renewal entropy estimator also becomes less  Figure 3: Simulated binary renewal process with characteristics similar to those of a bursty neuron. The ISI distribution is given by a mixture of discrete Gammas, 9 10 f (2,10) + 1 10 f (50,50) . The "rate" of the process is E(X 1 ) = 0.00398. The first plot shows the first 10000 simulated bits of X; the second plot shows the empirical ISI distribution; the last plot shows the autocorrelogram of X for values of the lag between 0 and 500. accurate. The situation is similar for the CTW. As shown in [12] [14], the CTW method also essentially approximates the empirical ISI distribution (although it does so in a different, more efficient way than the renewal entropy estimator), and for larger values of (α 2 , β 2 ) its estimates become less accurate, in a way analogous to the renewal entropy estimator. renewal plug-in (α 2 , β 2 ) entropy w = 20Ĥ nHn CTW (10,20) %  Table 6: Comparison of the ratios between the bias, the standard error and the square-root of the MSE over the true entropy rate. Five estimators are used on three different data sets generated by a renewal process whose ISI distribution P is a mixture of Gammas, P = µf (2,10) + (1 − µ)f (α 2 ,β 2 ) . The mixing parameter µ = 0.8 in the first two cases, and µ = 0.9 in the last case. All estimators are applied to data of the same total length n = 10 6 . Results are shown as percentages, that is, all ratios are multiplied by 100.

Summary and Concluding Remarks
A systematic and extensive comparison between several of the most commonly used and effective entropy estimators for binary time series was presented. Those were the plug-in or "maximum likelihood" estimator, four different LZ-based estimators, the CTW method, and the renewal entropy estimator.
Methodology. Three new entropy estimators were introduced; two new LZ-based estimators, and the renewal entropy estimator, which is tailored to data generated by a binary renewal process. A bootstrap procedure, similar to the one employed in [44], was described, for evaluating the standard error of the two sliding-window LZ estimatorsĤ n,k andH n,k . Also, for these two estimators, a practical rule of thumb was heuristically derived for selecting the values of the parameters n and k in practice.
Theory. It was shown (Theorem 2.1) that the two new LZ-based estimatorsH n,k andH n are universally consistent, that is, they converge to the entropy rate for every finite-valued, stationary and ergodic process. Unlike the corresponding estimatorsĤ n,k andĤ n of [22], no additional conditions are required. An effective method was derived (Proposition 3.1) for the accurate approximation of the entropy rate of a finite-state HMM with known distribution. Heuristic calculations were presented and approximate formulas derived for evaluating the bias and the standard error of each estimator.
Simulation. Several general conclusions can be drawn from the simulation experiments conducted.
(i) For all estimators considered, the main source of error is the bias.
(ii) Among all the different estimators, the CTW method is the most effective; it was repeatedly and consistently seen to provide the most accurate and reliable results.
(iii) No significant benefit is derived from using the finite-context-depth version of the CTW.
(iv) Among the four LZ-based estimators, the two most efficient ones are those with increasing window sizes,Ĥ n ,H n . No systematic trend was observed regarding which one of the two is more accurate.
(v) Interestingly (and somewhat surprisingly), in many of our experiments the performance of the LZ-based estimators was quite similar to that of the plug-in method.
(vi) The main drawback of the plug-in method is its computational inefficiency; with small word-lengths it fails to detect longer-range structure in the data, and with longer word-lengths the empirical distribution is severely undersampled, leading to large biases.
(vii) The renewal entropy estimator, which is only consistent for data generated by a renewal process, suffers a drawback similar (although perhaps less severe) to the plug-in.
In closing we note that much of the work reported here was done as part of the first author's Ph.D. thesis [12]. Several new estimators and results have appeared in the literature since then, perhaps most notably in [6]. There, a different entropy estimator is introduced based on the Burrows-Wheeler transform (BWT), it is shown to be consistent for all stationary and ergodic processes, and bounds on its convergence rate are obtained. Moreover, simulation experiments on binary data indicate that it significantly outperforms the plug-in estimator as well as a modified version of the LZ-based estimatorĤ n of [22]. In view of the present results, interesting further work would include a detailed comparison of the BWT estimator of [6] with the CTW method.