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.
Electrophysiological Classes of Cat Primary Visual Cortical Neurons In Vivo as Revealed by Quantitative Analyses (free article)
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
Bursts and recurrences of bursts in the spike trains of spontaneously active striate cortex neurons. (not free)
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
A nonparametric approach for detection of bursts in spike trains (not free)
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
Comparison of discharge variability in vitro and in vivo in cat visual cortex neurons (not free)
$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)