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)}
