Poisson Process Jumping between an Unknown Number of Rates: Application to Neural Spike Data
We introduce a model where the rate of an inhomogeneous Poisson process is modified by a Chinese restaurant process. Applying a MCMC sampler to this model allows us to do posterior Bayesian inference about the number of states in Poisson-like data. Our sampler is shown to get accurate results for synthetic data and we apply it to V1 neuron spike data to find discrete firing rate states depending on the orientation of a stimulus.
In reply to the author's feedback "Our description of uniformization might be misleading, because the problems we describe do not occur in all of its applications. For the SGCP model discussed in our paper, however, the uniformization really is over the rate, which is lambda in our model. There is no MJP in the SGCP model, because the rate is continuous." After rereading the relevant sections of the paper, I am sure that this is incorrect. From page 2 line 070 "in our model lambda(t) is a latent piecewise constant process. ... Let us assume we already have i segments [from the piecewise constant functino lambda(t)] and draw the next jump time from an exponential distribution with rate f". Therefore, the authors' process is an MJP with constant rate f, and state space given by the positive reals. The emission model of the MJP is a Poisson process with rate lambda(t). To use uniformization in this case, we must upper bound the rate of the transitions of the MJP --- that is done by chosing f. Surely, lambda(t) might be large during bursts, but f is still constant and so the number of auxiliary uniformization events does not depend on lambda(t). "For our neural data the maximum rate would have to be the spiking rate during the strongestbursts, but this would lead to a very large number of (later thinned out) event times to be sampled in the long periods between bursts, which slows down the algorithm severely." That is incorrect. The maximum rate used by uniformization is f (that's the transition rate of the latent piecewise linear process) which is constant and doesn't depend on the spike rate.---In the manuscript "Poisson Process Jumping between an Unknown Number of Rates: Application to Neural Spike Data" the authors provide a model for spike train data using an Markov jump process (MJP) modulating the rate of spikes. The state space of the MJP is given by a Dirichlet process (DP) supported on the rate of spikes. The DP is then marginalized, and Metropolis-Hastings updates are derived for the resulting CRP. The authors apply this model to simulated data and also bursting in V1 neurons of cats (data from NSF-funded CRCNS Data Sharing repository). Overall, I thought that the manuscript was good. I did however think that it could have made more contact with computational neuroscience literature. Also, the manuscript incorrectly describes uniformization and the implications for efficiency.The authors' approach is good: it applies modern machine learning to a very concrete problem. The manuscript is therefore topical for the Neural Information Processing Systems conference. However, I felt like this paper could have been more in touch with classical work on HMMs in computational neuroscience. For example, I felt that the authors should have extended their model to ensemble cases. HMM methods have long been used in computational neuroscience to model spike rates (see Abeles, M., & Gerstein, G. L. 1988. Detecting spatiotemporal firing patterns among simultaneously recorded single neurons. Journal of Neurophysiology, 60(909), 909–24. for example.) Rather than focusing on single neurons, each HMM state in Abeles & Gerstein 1988 is parameterized by a separate rate for each neuron in the ensemble. Almost every other work on HMMs in computational neuroscience also looks at ensemble cases.The same model presented . was used in Vinayak Rao's PhD thesis for modelling neural spikes (chapter 5.4.4 of MCMC for continuous-time discrete-state systems, University College London, 2012). In that work, a Dirichlet process prior was put on the spike rates of neurons in a grasshopper auditory cortex, and then the rates and transitions were learned using uniformization. Note that this work was not published outside of his thesis and also his thesis does not go into much detail about the paradigms used, so I would definitely not be opposed to anyone publishing a detailed account of it.The authors do mention uniformization as another approach to this sort of data (page 7 lines 351-365), but they state that it would be inefficient because uniformization requires an upper bound on the rate, and this upper bound would be large due to bursting. Unfortunately, this is a misunderstanding of the way uniformization works and it's not the case that uniformization would be inefficient, at least not for this reason. The upper bound required in uniformization is on the transition rate of the underlying MJP, not on the likelihood parameter (the spike rate). The author's model, viewed as an MJP, has transition rate f, which is a parameter they mention on page 2 line 78. Because f is constant over the process, an upper bound on the transition rate of the MJP is given by f, and this (not the spike rate) would be used as an upper bound in uniformization. I do agree that it is confusing because there are two rates here: (1) the transition rate of the latent piecewise linear function (which is f) and (2) the spike rate (which is lambda(t)), but the upper bound in uniformization cares only about (1) and does not care about (2).I found that the comment by the authors that the shifting, adding, removing and switching updates alone were "very unlikely to reduce the number of states" to be sort of strange. If those actions are balanced, then around a mode shouldn't they be just as likely to reduce the number of states as to increase the number of states? Is the issue that these updates alone won't mix and therefore won't get near a mode? Would be nice to make that more explicit.Another comment that I found sort of strange was "The amounts of both jumps and states seem to be captured well, but for a large number of distinct states the mean seems to underestimate the true value. This is not surprising, because the lambda parameters are drawn from the same base distribution. For a large number of states the probability that two states are very similar becomes high, which makes them indistinguishable without observing more data" (page 5 lines 256 to 260) which the authors wrote to explain why the emperical statistic of the posterior number of states in Figure 3 (left) falls below the y=x line in the qq-plot. While it is true that the model will be unable to distinguish between lambda values that are close together with a small amount of data, it is also true that due to sampling noise (MCMC variance) a sampler will also occasionally put data that is close together and well explained by one cluster into two (or more) separate cluster, and the rate at which this happens should sort of cancel out the rate of `incorrectly' clustering data from more than one (very similar) clusters into one cluster. Due to the main theorem behind MCMC, under the procedure: 1. sample parameters y from the prior, 2. sample data x|y, 3. sample Y|x through your MCMC scheme, the variable Y should be an unbiased estimator of y. So, unless I am missing something, it seems the mismatch in the qq-plot in Figure 3 (left) is still unexplained and could indicate failure to mix, or a bug.Minor notes:The citations within the text of this manuscript are not conforming to the NIPS standard. NIPS should cite references using [1], [2], [3] etc, and then have [1], [2], [3] preceeding the citations in the reference section. While this might seem pedantic, it should be changed for the camera ready copy should the paper be accepted, in order to allow the entire NIPS proceedings to maintain it's well deserved familiar and professional feeling.