# Estimating burstiness in a time series: Dependence on overall frequency?

I have a dataset of time series normalized for length. A working R example below. I am trying to estimate burstiness within the series, i.e. to compare series in the degree of clusteredness in these events. Here I just used a simple index of dispersion (variance of intervals divided by mean of intervals).

For this dataset it gives nice results that conform to intuitions, however they look too clean. With some exceptions, there seems to be a very strong dependency on the number of events within a series and their clusteredness. In fact we would expect the series with more events to have them less clustered. But I'm worried if this dependency we are seeing could be just simple mathematics, and nothing about the data.

Could this dependency be caused just by this simple way of calculation? Are there better ways to calculate clusteredness that would allow one to compare series with many events and series with few events - with both kinds having equal chance of being clustered and depending only on the distribution of events and not the amount?

This dependence for made up data is plotted below. The made up data is not very diagnostic but shows a similar pattern to overall data.

#make grid
grid = rbind(rep(0, 1000),rep(0, 1000),rep(0, 1000))
#add 3 event series on grid
grid[1,c(0, 35, 71, 179, 181, 229, 244, 249, 269, 270, 279, 281, 297, 311, 322, 335, 342, 415, 420, 421, 433, 435, 436, 437, 459, 462, 469, 549, 551, 552, 556, 559, 583, 804, 813, 860, 861, 864, 867, 868, 869, 1000)]<-1
grid[2,c(0, 35, 297, 311, 322, 335,  437, 459, 462, 469, 549, 551, 552, 804, 813, 860, 861, 864, 867, 868, 869, 1000)]<-1
grid[3,c(0, 229, 244, 249, 269, 270, 279, 552, 556, 559, 583, 1000)]<-1
#output variance to mean ratios and number of events for each series
VMRs <- rep(NA,length(grid[,1]))
events <- rep(NA,length(grid[,1]))
for(j in 1:length(grid[,1])){
tlength <- length(which(grid[j,]==1))
intervals <- rep(NA,tlength-1)
for(i in 2:tlength){
intervals[i-1] = (which(grid[j,]==1)[i]-which(grid[j,]==1)[i-1])/1000
}
VMRs[j]=var(intervals)/mean(intervals)
events[j]=tlength
}
#plot the results
plot(VMRs,events)


EDIT: I simulated some data supposing a uniform chance of each point becoming an event, and there definitely is a dependency with frequency of events in this simple equation. This may be related to the fact that the events take place in a limited time? Yes, more chances to interrupt the intervals probably gives a more even distribution, but there ought to be a measure that takes this into account. I have trouble comprehending why exactly the number of events has an influence. Any advice will be welcome!

VMRs <- matrix(NA,100,100)#rep(NA,100)
for (events in 1:100){
for (j in 1:100){
grid = rep(0, 1000)
grid[c(floor(runif(events)*1000),1,1000)]<-1
intervals <- rep(NA,length(which(grid==1))-1)
for(i in 2:length(which(grid==1))){
intervals[i-1] = (which(grid==1)[i]-which(grid==1)[i-1])/1000
}
VMRs[j,events]=var(intervals)/mean(intervals)
}
}
colMeans(VMRs)

plot(colMeans(VMRs))