Spike train entropy-rate estimation using hierarchical Dirichlet process priors
Entropy rate quantifies the amount of disorder in a stochastic process. For spiking neurons, the entropy rate places an upper bound on the rate at which the spike train can convey stimulus information, and a large literature has focused on the problem of estimating entropy rate from spike train data. Here we present Bayes least squares and empirical Bayesian entropy rate estimators for binary spike trains using hierarchical Dirichlet process (HDP) priors. Our estimator leverages the fact that the entropy rate of an ergodic Markov Chain with known transition probabilities can be calculated analytically, and many stochastic processes that are non-Markovian can still be well approximated by Markov processes of sufficient depth. Choosing an appropriate depth of Markov model presents challenges due to possibly long time dependencies and short data sequences: a deeper model can better account for long time dependencies, but is more difficult to infer from limited data. Our approach mitigates this difficulty by using a hierarchical prior to share statistical power across Markov chains of different depths. We present both a fully Bayesian and empirical Bayes entropy rate estimator based on this model, and demonstrate their performance on simulated and real neural spike train data.
The paper introduces a way to estimate the entropy 
      rates (and transition probabilities) of a Markov Model from 
      observations, by improving the empirical bin count through a 
      tree-based Bayesian prior which takes into account the progress of the 
      probabilities through time and assigns more similar future 
      probabilities to observations which differ further in the past. 
      The paper has some interesting ideas, and the sequential 
      construction of the Dirichlet priors and the extraction of the 
      estimates is elegant. Certainly, the results look good - but in 
      Fig. 3(c), the superiority of convergence of the HDP estimators is 
      not really clear compared to the others. The text in lines 407-409 
      seems to insinuate that the performance of the HDP estimators is 
      somehow special whereas it is quite on coherence with the other 
      methods. The transition probability, of course, is nicely 
      regularized, but stronger data is required to have an idea how good 
      this regularization is. Fig. 2(a), for instance shows a good match 
      of the high-probability transitions, but the low-probability 
      transitions seem not as well matched. This is slightly surprising, 
      as there are only two possibilities for the transition, so it is not 
      clear why the high transition probability values should be 
      regularized to a higher degree than the low transition probability 
      values. nevertheless, one may wish some more reference to 
      other reconstruction methods, such as CSSR by Shalizi et al, 
      suffix-tree methods (Y. Singer et al.), observable operator methods 
      (Jaeger). There is a lot of emphasis on the weaknesses of pure 
      entropy estimation in the paper, but clearly also the structure of 
      the reconstruction mechanism hovers in the background which deserves 
      mention as your construction follows the spirit of the epsilon 
      machine reconstruction (Shalizi & Crutchfield) and may find 
      relevant generalizations in that direction. Detailed comments: 
      - line 124: what's e in log(e)? - line 128: second H_NSB 
      has wrong subscript - line 177: it is not clear what the bold 
      characters mean - throughout the paper, until quite late, I was 
      confused whether the new symbol appears at the *right* or the *left* 
      end of the string, and my confusion was caused mainly by this 
      illustration. Why does 010 transitionf to 101? If the new symbol 
      is on the right, where does the 0 at the left end disappear to? Do 
      we have a 2-time step memory system only here? - line 227: 
      similar issue - "g_0011 should be similar to g_011" - which one is the 
      latest symbol. Here it seems it may be the first one. Later in the 
      text you seem to refer to the latest symbol at the right. Please 
      clarify. - line 249: "whose sample paths are each probability 
      distributions" - this sentence is unclear, information seems 
      missing; what condition do the probability distributions 
      associated with the paths have to fulfil? - line 264: there is 
      information missing here - it seems that when you start with a fresh 
      distribution (p_0), you still seem to assume existence of N 
      observations (the length of your data string). How come you do not 
      start with the empty string straight away and build up your Dirichlet 
      prior as you go through the first N entries? - line 269: I do 
      not understand the structure of the posterior predictive distribution. 
      You state the set g_{|k|}, but I do not see the connection between the 
      set and the required structure of the distribution. Reformulate. Also 
      the rest of the paragraph below line 272 is unclear. Please fix that, 
      maybe it is a consequence that I did not understand your construction 
      of g_pr. - line 294: I have the impression that the line contains 
      an error: the Beta function for the g_{0s} and the g_{1s} terms 
      probably need an \alpha_{|s|+1} instead of \alpha_{|s|} - line 
      301: "each \alpha as fixed hyperparameters" - elide the "s"