# Inner correlation of occurrences (burstiness?) in R

I want to measure the inner correlation of the occurrences of events i.e. I want to distinguish between the two (drawn) and say "in the second sample events occur more conglomerate compared to the first":

Isn't this different from burstiness (i.e. easier to compute as it doesn't invole traffic-values)?

I have drawn a possible representation of the expected result (for S2) in the image below the samples (just as it is in my mind now), but don't hesitate to suggest different proposals.

I looked out and tried different build-in functions and googled/ searched CrossValidated for R and burstiness or "inner correlation", but I didn't make progress. Maybe I am looking for the wrong search terms. If I wanted to really measure burstiness, it seemed to me (by reading papers about burstiness) that there would be no common known measure (but many different one). I we would pick one it would be cool to justify why choosing this specific one.

This is how the time series is represented in R, currently: (written by dput in file)

c(3.256861, 3.377142, 3.941173, 4.304236, 4.485358, 4.606512,
4.707296, 5.473004, 5.714746, 5.815394, 5.835405, 5.936067, 5.957008,
6.964611, 7.045158, 7.065171, 7.165824, 7.669618, 8.17324, 8.273692,
9.503988, 9.604991, 9.624853, 9.725522, 10.237766, 10.954529,
11.378399, 12.687714, 13.291919, 13.41258, 13.67527, 14.380529,
14.743638, 15.247138, 15.851832, 15.952875, 15.972497, 16.456259,
16.476052, 17.201506, 17.463708, 18.068535, 18.309645, 18.390292,
18.410299, 18.430323, 18.531736, 18.652921, 18.793662, 19.297076,
19.639692, 19.760698, 20.768096, 20.868441, 20.990499, 21.494412,
21.856368, 22.199341, 22.219143, 22.440472, 22.481118, 23.327013,
23.447678, 23.811188, 23.843, 24.113302)


This is a widely studied problem in neurosciences, where you need to determine the "burstiness" of action potentials of a neuron. The methods, however, can be obviously applied to any series of events.

Most of them rely on the analysis of the intervals between two following events: in the case of action potentials these are generally called inter-spike interval (ISI), but we can call them inter-event intervals (IEI) to generalize.

We can define them as

$IEI = t_n - t_{n-1} \quad\quad n=2,3,4,...,N$

Where $t_n$ is the time of event $n$ and $N$ is the total number of events.

I will list some of the approaches that have been used. Mind, however, that this list is far from exaustive.

The easiest visual thing to do is starting to plot an histogram of the IEIs or, even better, an histogram of $log_{10}(IEI)$. In case of high "burstiness" the histogram will have a clear bimodal distribution, with lots of short intervals between events and some longer ones (the pauses between bursts)

If you have a fairly good number of series of event you can also use a clustering algorithm to divide them in groups (regular, slow bursting, fast bursting etc.). This approach was taken, for instance, in this paper by Nowak et al. where several parameters of the distribution (mean, median, skewness, kurtosis, IQI etc.) are taken as classifiers for hierarchical clustering.

Another classic approach is known as the "Poisson surprise method" and was described in 1985 by Charles Legéndy and Michael Salcman in their paper

The idea of the method is that:

The measure used here is an evaluation of how improbable it is that the burst is a chance occurrence and is computed, for any given burst that contains n spikes in a time interval T, as

$s = - log P$

where P is the probability that, in a random (Poisson) spike train having the same average spike rate Y as the spike train studied, a given time interval of length T contains y2 or more spikes.

I can provide R code for this if needed

An "updated" version of the Poisson-surprise method, which was developed to solve certain issues with that method is the rank-surprise method described in 2007 by Boris Gourévitch and Jos Eggermont in their paper

which uses a non parametric approach to define bursts.

We propose to use a more exhaustive search of the maximum of the surprise statistic using the following algorithm dubbed ESM (exhaustive surprise maximization): preliminary to the algorithm, we fix the largest ISI value acceptable in a burst (limit) and a level −log(α) of minimum significance for the surprise statistic. We, then identify a first sequence of ISIs whose values are below limit. From this sequence, we perform an exhaustive search of the highest surprise statistic over all possible continuous subsequences of ISI.If the final surprise statistic is above −log(α), the associated subsequence is labeled as a burst. Another burst is then searched among the remaining continuous ISI subsequences, obeying the same criterion. The process is repeated until one remaining continuous ISI subsequence is able to provide a significant RS statistic. When the process stops, we proceed to the next sequence of ISIs whose values are below limit and so on.

The authors provide pseudo-code and Matlab code for the algorithm

Other approaches rely on the variability of the distribution.

In particular, one can use the coefficient of variation $C_V$, classically defined as

$C_V = \frac{\sigma_{IEI}}{\langle{IEI}\rangle}$

The higher $C_V$ the burstier the events' distribution

$C_V$, however, is a fairly rough index, so a finer version of it was proposed, called $C_{V2}$, by Gary Holt and colleagues in their paper

$C_{V2} = \frac{2*|{IEI}_{n+1}-{IEI}_n|}{{IEI}_{n+1}+{IEI}_n}$

Finally, another approach, proposed by Shigeru Shinomoto and colleagues in 2003 is the local variation coefficient $L_v$ which is defined as

$L_v = \frac{1}{n-1} \sum_{i=1}^{n-1}\frac{3(T_i-T_{i+1})^2}{(T_i+T_{i+1})^2}$

in their paper

Differences in Spiking Patterns Among Cortical Neurons

Also, two classical must-reads:

Neuronal spike trains and stochastic point processes. I. The single spike train

Neuronal spike trains and stochastic point processes. II. Simultaneous spike trains (both free, the second one probably is not too interesting for you, but it's still a good read)

• Thank you nico, your detailed answer helped a lot. I am not sure if the "Poisson surprise method" is that intuitive, so I will propaply go for the $C_{V2}$, by Gary Holt et.al. which I don't need to describe to long (and I cannot see a straight advantage) and definitely make sense to me. Furthermore the hint towards a cluster analysis is awesome, as I really have a whole lot of series. – danielberger Jun 19 '11 at 21:33
• @danielberger: you're welcome! I forgot to say that indices like Cv2 and Lv are "burstiness indices" that tell you how bursty the time series is. Methods like Poisson/rank surprise are instead "burst identification methods", so you end up with a list of bursts. For instance the method will tell that event 1 is alone, then you have a burst between event 2 and 10 and then 3 separate events, etc etc. Finally, I find the Lv index to be more accurate than Cv2 in some situations but it really depends on your data. If the data is long enough you can also try to calculate it on a sliding window. – nico Jun 20 '11 at 6:20
• @danielberger: also, if you don't have access to any of those papers I can send them to you in PDF. – nico Jun 20 '11 at 6:23