Non-parametric (smoothed) estimate of current rate

Question

I am looking at time-series event data and need to visualize how the arrival rate $\lambda$ changes over time. I do not want to assume any underlying distribution that the data comes from (it is certainly not a Poisson process as the rate changes, and neither the Weibull nor the log-logistic distribution seem to describe what's going on). What I want is something between the following two extremes:

1. Assume that the rate does not change with the time since the previous event (Poisson), and simply calculate the overall rate of events per time. Problem: just a single number, i.e., no ability for the rate to change over time.
2. Calculate the empirical instantaneous rate between all sets of neighboring points (as in $\frac{1}{t_i - t_{i-1}}$) and plot this over time. Problem: too jumpy, not very informative.

I want something in between these two extremes, in the sense that a kernel density estimate is "in between" calculating the mean of observations, vs. looking at a rugplot showing the individual values of the data. Ultimately, I want to see a smoothed curve that shows a locally-weighted average of the rate of arrivals over time.

Two classes of approach come to mind: some kind of kernel density estimate of the values from extreme #2, or some use of an (inverse) exponential function of the previous interarrival times to provide an exponentially-weighted average of the rate.

Is there a right way to do this? If so, is there also a right way to select an appropriate bandwidth/weighting value (which I'm assuming will be necessary, whatever the solution is)?

PS This is not a failure analysis, but the concept of hazard rates seems somewhat relevant. However, I'm not finding much on how to get the hazard rate function without assuming a distribution.

Answer 1

This is a common problem in neuroscience, and the same methods can be applied to your problem. I'll describe the neuroscience problem briefly, because this will help to interpret the papers below. The methods I'll describe are covered in these papers.

In order to transmit information, neurons fire action potentials (a.k.a. spikes), which are sharp, transient increases in voltage across the cell membrane. Experiments often measure spikes using electrodes, and the resulting time series is called a spike train. A spike train is characterized by the times the spikes occurred, and can be modeled as a point process. We're often interested in the instantaneous firing rate (i.e. the expected number of spikes fired per unit time), which typically fluctuates over time. This must be estimated from the recorded spike times. Spikes are analogous to your events, and instantaneous firing rates are analogous to your time-varying event rate.

The simplest method for estimating instantaneous firing rates is called kernel smoothing. Here, firing rates are obtained by convolving the spike train with a kernel/window function. Suppose spikes occur at times $S = \{s_1, \dots, s_n\}$. Let $k_\sigma(\cdot)$ be a kernel function whose width is controlled by the bandwidth parameter $\sigma$. Then the estimated firing rate at time $t$ is:

$$f(t) = \sum_{i=1}^n k_\sigma(t - s_i)$$

The principle here is the same as in kernel density estimators (KDEs). As with KDEs, the firing rate estimate does not depend strongly on the shape of the kernel, but does depend strongly on the bandwidth. Larger bandwidths give smoother estimates that are less sensitive to noise, but cannot capture fast changes. Smaller bandwidths give less smooth estimates that better capture fast changes, but are more sensitive to noise. The bandwidth is often chosen in an ad hoc manner, or using prior knowledge/assumptions about the underlying signal. Methods have also been developed for optimal bandwidth selection, by modeling the spike train as an inhomogeneous Poisson process and treating the kernel smoothed firing rate as an estimate of the underlying rate function. Other methods allow the bandwidth itself to vary over time.

More complex Bayesian estimators have also been developed. Here, a prior is defined over the underlying firing rate function, e.g. using Gaussian processes or splines. The prior encodes assumptions about how the firing rates vary over time, such as smoothness. A likelihood function specifies how spikes are generated, given the underlying firing rate. This takes the form of a non-stationary point process (e.g. an inhomogeneous Poisson process or Gamma interval process). Bayesian inference is then used to obtain a posterior distribution over the firing rate, given the observed spike times.

References

Cunningham, J. P., Gilja, V., Ryu, S. I., & Shenoy, K. V. (2009). Methods for estimating neural firing rates, and their application to brain–machine interfaces. Neural Networks, 22(9), 1235-1246.

Tomar, R. (2019). Methods of firing rate estimation. Biosystems, 183, 103980.