# Entropy estimation for a symbol sequence

I am looking for an R-implementation of the Lempel-Ziv data compression algorithm, to estimate the source entropy of a time-series consisting of a sequence of symbols.

Rather than simply measuring the entropy of the (time-aggregated) symbol-frequency distribution, the Lempel-Ziv algorithm can be used to estimate the entropy of the temporal ordering within the sequence itself.

The sequences I am studying are derived from animal trajectories; each $x,y$ location is assigned an arbitrary label (A,B,...,Z). Represented in this way, the trajectory of a moving animal looks like B,G,G,S,Y,G,Y,H,D,S,A,B,B,G,G,G,H etc. Hence, I am interested in measuring the $\textit{predictability}$ of the sequence.

The procedure for estimating the source entropy is (briefly) described in the Supplementary Information of Song et al (2010), and they describe the Lempel-Ziv estimator as:

$$S^{est} = \left( \frac {1}{n} \sum_{i=1} \Lambda_i \right)^{-1} ln~n$$ where $\Lambda_i$ is the length of the shortest substring starting at position $i$ which $\textit{doesn't}$ previously appear from position 1 to $i$-1.

So, does anyone have any suggestions for how to compute the above in R?

Incidentally, I should mention that I have multiple trajectories between which I would like to make entropy comparisons, however, some of the trajectories contain gaps that correspond to times when the animal was not observed.

REFERENCE:

Song, Qu, Blumm & Barabasi (2010) Limits of Predictability in Human Mobility. $\textit{Science}$. 19. http://www.sciencemag.org/content/327/5968/1018.short

UPDATE: The function below is what I arrived at. However, rather than chugging through all the possible sub-sequences after the i'th index, and then finding the shortest novel sub-sequence, it would be faster to use a while loop. Similarly I am updating the dictionary of newly-acquired unique substrings using rbind, which is not efficient. OTOH the sequences I'm using are short (<300), so this isn't a problem, but there is certainly room for improvement.

PREDICTABILITY <- function(Sequence, Max_String)
{

output     <- matrix(NA,nrow=length(Sequence), ncol=1)
dictionary <- NULL
## cycle through each entry, i, in the Sequence
for (i in 1:(length(Sequence)-Max_String) )
{
## Compile list of increasingly-long substrings, starting at position i; i+0,i+1,i+2,...,i+Max_String
codons <- matrix(NA, nrow=Max_String, ncol=1)
for (STRL in 0:Max_String)
{
codons[STRL,] <- paste(Sequence [i:(i+STRL-1)], collapse="")
}
## Find which of these codons have NOT been seen before
new <- codons[!codons %in% dictionary]
## check for no new codons
ifelse ((length(new)>0),
record <- min(nchar(new)),    ## if we have new codons, find the shortest among them
record <- NA )                ## if none are new (because we aren't searching far enough ahead), assign NA...
## find the shortest of these unseen codons
output[i,] <- record

## Finally, add the unseen codons to the dictionary
dictionary <- c(dictionary, new)
}##i
## Calculate source entropy (i.e. predictability) from formula in Song et al (2010)
n <- length(output[!is.na(output)])
S <- (1/mean(output, na.rm=TRUE)) * log(n)  ## Source entropy ?natural log means the units are nats?
return(S)}


I can't quite understand the PREDICTABILITY function you introduced.

Here's some of experiments I did.

PREDICTABILITY(rep(c(0,1,2,3),250),10)

[1] 1.386294

PREDICTABILITY(floor(runif(1000,min=0,max=4)),10)

[1] 1.339492

PREDICTABILITY(rep(0,1000),10)

[1] 0

From rep(0,1000), I figure it's the measure of complexity of a sequence. But by comparison of rep(c(0,1,2,3),250) and floor(runif(1000,min=0,max=4)), I am quite confused. If it's a measure of complexity of a sequence, the measure should be higher for random sequence(runif) then repetition of 0,1,2,3.

Here are other examples.

PREDICTABILITY(rep(0:1,500),10)

[1] 0.6931472

PREDICTABILITY(floor(runif(1000,min=0,max=2)),10)

[1] 0.7573678

PREDICTABILITY(rep(c(0,1,1,0),250),10)

[1] 0.9241962