Abstract

Kernel smoother and a time-histogram are classical tools for estimating an instantaneous rate of spike occurrences. We recently established a method for selecting the bin width of the time-histogram, based on the principle of minimizing the mean integrated square error (MISE) between the estimated rate and unknown underlying rate. Here we apply the same optimization principle to the kernel density estimation in selecting the width or "bandwidth" of the kernel, and further extend the algorithm to allow a variable bandwidth, in conformity with data. The variable kernel has the potential to accurately grasp non-stationary phenomena, such as abrupt changes in the firing rate, which we often encounter in neuroscience. In order to avoid possible overfitting that may take place due to excessive freedom, we introduced a stiffness constant for bandwidth variability. Our method automatically adjusts the stiffness constant, thereby adapting to the entire set of spike data. It is revealed that the classical kernel smoother may exhibit goodness-of-fit comparable to, or even better than, that of modern sophisticated rate estimation methods, provided that the bandwidth is selected properly for a given set of spike data, according to the optimization methods presented here.

Fixed kernel density estimation. (a) The underlying spike rate λt of the Poisson point process. (b) 20 spike sequences sampled from the underlying rate, and (c) Kernel rate estimates made with three types of bandwidth: too small, optimal, and too large. The gray area indicates the underlying spike rate. (d) The cost function for bandwidth w. Solid line is the estimated cost function, Eq. (), computed from the spike data; The dashed line is the exact cost function, Eq. (), directly computed by using the known underlying rate

Variable kernel density estimation. (a) Kernel rate estimates. The solid and dashed lines are rate estimates made by the variable and fixed kernel methods for the spike data of Fig. (b). The gray area is the underlying rate. (b) Optimized bandwidths. The solid line is the variable bandwidth determined with the optimized stiffness constant γ ∗ = 0.8, selected by Algorithm 2; the dashed line is the fixed bandwidth selected by Algorithm 1. (b) The cost function for bandwidth stiffness constant. The solid line is the cost function for the bandwidth stiffness constant γ, Eq. (), estimated from the spike data; the dashed line is the cost function computed from the known underlying rate

Fitting performances of the six rate estimation methods, histogram, fixed kernel, variable kernel, Abramson’s adaptive kernel, Locfit, and Bayesian adaptive regression splines (BARS). (a) Two rate profiles (2 [s]) used in generating spikes (gray area), and the estimated rates using six different methods. The raster plot in each panel is sample spike data (n = 10, superimposed). (b) Comparison of the six rate estimation methods in their goodness-of-fit, based on the integrated squared error (ISE) between the underlying and estimated rate. The abscissa and the ordinate are the ISEs of each method applied to sinusoidal and sawtooth underlying rates (10 [s]). The mean and standard deviation of performance evaluated using 20 data sets are plotted for each method

Application of the fixed and variable kernel methods to spike data of an MT neuron (j024 with coherence 51.2% in nsa2004.1 (Britten et al. )). (a–c): Analyses of n = 1, 10, and 30 spike trains; (top) Spike rates [spikes/s] estimated with the fixed and variable kernel methods are represented by the gray area and solid line, respectively; (middle) optimized fixed and variable bandwidths [s] are represented by dashed and solid lines, respectively; (bottom) A raster plot of the spike data used in the estimation. (d) Comparison of the two kernel methods. Bars represent the difference between the cross-validated cost functions, Eq. (), of the fixed and variable kernel methods (fix less variable). The positive difference indicates superior fitting performance of the variable kernel method. The cross-validated cost function is obtained as follows. Whole spike sequences () are divided into blocks, each composed of n ( = 1, 5, 20, 30) spike sequences. A bandwidth was selected using spike data of a training block. The cross-validated cost functions, Eq. (), for the selected bandwidth are computed using the leftover test blocks, and their average is computed. The cost function is repeatedly obtained -times by changing the training block. The mean and standard deviation, computed from samples, are displayed