5
$\begingroup$

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)}
$\endgroup$
1
$\begingroup$

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.