J Comput Neurosci (2010) 29:171–182
DOI 10.1007/s10827-009-0180-4

Kernel bandwidth optimization in spike rate estimation
Hideaki Shimazaki · Shigeru Shinomoto

Received: 22 December 2008 / Revised: 20 May 2009 / Accepted: 23 July 2009 / Published online: 5 August 2009
© The Author(s) 2009. This article is published with open access at Springerlink.com

Abstract Kernel smoother and a time-histogram areclassical tools for estimating an instantaneous rate of
spike occurrences. We recently established a method
for selecting the bin width of the time-histogram, basedon the principle of minimizing the mean integrated
square error (MISE) between the estimated rate andunknown underlying rate. Here we apply the same
optimization principle to the kernel density estimationin selecting the width or “bandwidth” of the kernel,
and further extend the algorithm to allow a variablebandwidth, 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 toavoid possible overfitting that may take place due toexcessive freedom, we introduced a stiffness constantfor bandwidth variability. Our method automaticallyadjusts the stiffness constant, thereby adapting to theentire set of spike data. It is revealed that the classicalkernel smoother may exhibit goodness-of-fit comparable to, or even better than, that of modern sophisticatedrate estimation methods, provided that the bandwidthis selected properly for a given set of spike data, according to the optimization methods presented here.

Action Editor: Peter Latham
H. Shimazaki (B)Grün Unit, RIKEN Brain Science Institute,
Saitama 351-0198, Japane-mail: shimazaki@brain.riken.jp
S. Shinomoto
Department of Physics, Kyoto University,
Kyoto 606-8502, Japane-mail: shinomoto@scphys.kyoto-u.ac.jp

Keywords Kernel density estimation ·Bandwidth optimization · Mean integrated
squared error

1 Introduction
Neurophysiologists often investigate responses of asingle neuron to a stimulus presented to an animal
by using the discharge rate of action potentials, or
spikes (Adrian 1928; Gerstein and Kiang 1960; Abeles
1982). One classical method for estimating spike rate is
the kernel density estimation (Parzen 1962; Rosenblatt
1956; Sanderson 1980; Richmond et al. 1990; Nawrot
et al. 1999). In this method, a spike sequence is convoluted with a kernel function, such as a Gauss density
function, to obtain a smooth estimate of the firingrate. The estimated rate is sometimes referred to as a
spike density function. This nonparametric method isleft with a free parameter for kernel bandwidth that
determines the goodness-of-fit of the density estimate
to the unknown rate underlying data. Although theories have been suggested for selecting the bandwidth,
cross-validating with the data (Rudemo 1982; Bowman
1984; Silverman 1986; Scott and Terrell 1987; Scott
1992; Jones et al. 1996; Loader 1999a, b), individual
researchers have mostly chosen bandwidth arbitrarily.
This is partly because the theories have not spreadto the neurophysiological society, and partly due toinappropriate basic assumptions of the theories themselves. Most optimization methods assume a stationaryrate fluctuation, while the neuronal firing rate often
exhibits abrupt changes, to which neurophysiologists, inparticular, pay attention. A fixed bandwidth, optimized

172

using a stationary assumption, is too wide to extract thedetails of sharp activation, while in the silent period, thefixed bandwidth would be too narrow and may cause
spurious undulation in the estimated rate. It is therefore
desirable to allow a variable bandwidth, in conformity
with data.
The idea of optimizing bandwidth at every instant was proposed by Loftsgaarden and Quesenberry(1965). However, in contrast to the progress in methods
that vary bandwidths at sample points only (Abramson
1982; Breiman et al. 1977; Sain and Scott 1996; Sain
2002; Brewer 2004), the local optimization of bandwidth at every instant turned out to be difficult because of its excessive freedom (Scott 1992; Devroyeand Lugosi 2000; Sain and Scott 2002). In earlier
studies, Hall & Schucany used the cross-validation of
Rudemo and Bowman, within local intervals (Hall and
Schucany 1989), yet the interval length was left free.
Fan et al. applied cross-validation to locally optimized
bandwidth (Fan et al. 1996), and yet the smoothness ofthe variable bandwidth was chosen manually.
In this study, we first revisit the fixed kernel method,and derive a simple formula to select the bandwidth of
the kernel density estimation, similar to the previous
method for selecting the bin width of a peristimulus
time histogram (See Shimazaki and Shinomoto 2007).
Next, we introduce the variable bandwidth into the kernel method and derive an algorithm for determining the
bandwidth locally in time. The method automaticallyadjusts the flexibility, or the stiffness, of the variablebandwidth. The performance of our fixed and variable
kernel methods are compared with established density estimation methods, in terms of the goodness-offit to underlying rates that vary either continuously or
discontinuously. We also apply our kernel methods tothe biological data, and examine their ability by crossvalidating with data.
Though our methods are based on the classical kernel method, their performances are comparable to various sophisticated rate estimation methods. Because of
the classics, they are rather convenient for users: themethods simply suggest bandwidth for the standard
kernel density estimation.

J Comput Neurosci (2010) 29:171–182

The recorded spike trains are aligned at the onset ofstimuli, and superimposed to form a raw density, as
1
xt =δ (t − ti ),
n i=1
N

(1)

where n is the number of repeated trials. Here, each
spike is regarded as a point event that occurs at an
instant of time ti (i = 1, 2, · · · , N) and is represented
by the Dirac delta function δ(t). The kernel density
estimate is obtained by convoluting a kernel k(s) to theraw density xt ,(2)
λ̂t = xt−s k(s) ds.

Throughout this study, thethat does not ∞ integral
specify bounds refers to −∞ . The kernel function satisfies the normalization
condition, k(s) ds = 1, a zerofirst moment,
sk(s)
ds
=
0, and has a finite bandwidth, 22
w = s k(s) ds < ∞. A frequently used kernel is theGauss density function,1s2
kw (s) = √(3)
exp − 2 ,2w
2πw
where the bandwidth w is specified as a subscript. In the
body of this study, we develop optimization methods
that apply generally to any kernel function, and derivea specific algorithm for the Gauss density function in
the Appendix.
2.2 Mean integrated squared error optimization
principle
Assuming that spikes are sampled from a stochastic
process, we consider optimizing the estimate λ̂t to beclosest to the unknown underlying rate λt . Among
several plausible optimizing principles, such as the
Kullback-Leibler divergence or the Hellinger distance,
we adopt, here, the mean integrated squared error(MISE) for measuring the goodness-of-fit of an estimate to the unknown underlying rate, as bE(λ̂t − λt )2 dt,(4)
MISE =a2 Methods
2.1 Kernel smoothing
In neurophysiological experiments, neuronal responseis examined by repeatedly applying identical stimuli.

where E refers to the expectation with respect to the
spike generation process under a given inhomogeneous
rate λt . It follows, by definition, that Ext = λt .
In deriving optimization methods, we assume the
Poisson nature, so that spikes are randomly sampled at
a given rate λt . Spikes recorded from a single neuron
correlate in each sequence (Shinomoto et al. 2003, 2005,2009). In the limit of a large number of spike trains,

J Comput Neurosci (2010) 29:171–182

173

By noting that λt = Ext , the integrand of the second
term in Eq. (5) is given asExt Eλ̂t = E xt λ̂t − E (xt − Ext ) λ̂t − Eλ̂t ,(6)from a general decomposition of covariance of two
random variables. Using Eq. (2), the covariance (thesecond term of Eq. (6)) is obtained as
E (xt − Ext ) λ̂t − Eλ̂t= kw (t − s) E [(xt − Ext ) (xs − Exs )] ds
==however, mixed spikes are statistically independent andthe superimposed sequence can be approximated as asingle inhomogeneous Poisson point process (Cox 1962;
Snyder 1975; Daley and Vere-Jones 1988; Kass et al.
2005).

kw (t − s) δ (t − s)1
kw (0)Ext .
n

Cn (w) = MISE −
=

 b
a

 b
a

a

−2

 b

(5)

a

Rudemo and Bowman suggested the leave-one-out
cross-validation to estimate the second term of Eq. (5)(Rudemo 1982; Bowman 1984). Here, we directly estimate the second term with the Poisson assumption (Seealso Shimazaki and Shinomoto 2007).



1
E xt λ̂t − kw (0)Ext
n


dt.

(8)

Equation (8) is composed of observable variables only.
Hence, from sample sequences, the cost function isestimated as

 b
 b
1
2
Ĉn (w) =xt λ̂t − kw (0)xt dt.
(9)λ̂t dt − 2naa
In terms of a kernel function, the cost function is written as
1 Ĉn (w) = 2
ψw ti , t j
n i, j
⎫
⎧
⎬
2 ⎨
− 2
kw ti − t j − kw (0)N
⎭
n ⎩
i, j

=

λt Eλ̂t dt.

 ba

λ2t dtEλ̂2t dt − 2(7)Here, to obtain the second equality, we used the assumption of the Poisson point process (independent
spikes).
Using Eqs. (6) and (7), Eq. (5) becomes b
Cn (w) =Eλ̂2t dt

2.3 Selection of the fixed bandwidth
Given a kernel function such as Eq. (3), the density
function Eq. (2) is uniquely determined for a raw density Eq. (1) of spikes obtained from an experiment.
A bandwidth w of the kernel may alter the density
estimate, and it can accordingly affect the goodness-offit of the density function λ̂t to the unknown underlying
rate λt . In this subsection, we consider applying a kernelof a fixed bandwidth w, and develop a method forselecting w that minimizes the MISE, Eq. (4).
The integrand of the MISE is decomposed into threeparts: Eλ̂2t − 2λt Eλ̂t + λ2t . Since the last component
does not depend on the choice of a kernel, we subtract itfrom the MISE, then define a cost function as a functionof the bandwidth w:

1
Exs ds
n

1 
n2

i, j

ψw ti , t j −2 
kw ti − t j ,n2 i= j

(10)

where
ψw ti , t j = bakw (t−ti ) kw t−t j dt.

(11)

The minimizer of the cost function, Eq. (10), is an
estimate of the optimal bandwidth, which is denoted byw ∗ . The method for selecting a fixed kernel bandwidthis summarized in Algorithm 1. A particular algorithm

174

J Comput Neurosci (2010) 29:171–182

with data. The spike rate estimated with the variablebandwidth wt is given by(12)
λ̂t = xt−s kwt (s) ds.
Here we select the variable bandwidth wt as a fixedbandwidth optimized in a local interval. In this approach, the interval length for the local optimization
regulates the shape of the function wt , therefore, it
subsequently determines the goodness-of-fit of the estimated rate to the underlying rate. We provide a method
for obtaining the variable bandwidth wt that minimizes
the MISE by optimizing the local interval length.
To select an interval length for local optimization, we
introduce the local MISE criterion at time t as


2
u−t
localMISE = E λ̂u − λu ρW
du,
(13)

where λ̂u = xu−s kw (s) ds is an estimated rate with au−t
fixed bandwidth w. Here, a weight function ρW
localizes the integration of the squared error in a characteristic interval W centered at time t. An example of theweight function is once again the Gauss density function. See the Appendix for the specific algorithm for
the Gauss weight function. As in Eq. (5), we introduce
the local cost function at time t by subtracting the term
irrelevant for the choice of w as

Cnt (w, W) = localMISE − λ2u ρ u−t
(14)
W du.
The optimal fixed bandwidth w ∗ is obtained as a minimizer of the estimated cost function:
1  tti , t jĈnt (w, W) = 2
ψn i, j w,W
2 
ti −t
− 2
kw ti − t j ρW,
n i= j
where
t
ψw,W
ti , t j =developed for the Gauss density function is given in the
Appendix.
2.4 Selection of the variable bandwidth
The method described in Section 2.3 aims to selecta single bandwidth that optimizes the goodness-of-fitof the rate estimate for an entire observation interval[a, b ]. For a non-stationary case, in which the degreeof rate fluctuation greatly varies in time, the rate estimation may be improved by using a kernel function
whose bandwidth is adaptively selected in conformity

(15)


u−t
kw (u − ti ) kw u − t j ρW
du.

(16)

The derivation follows the same steps as in the previous section. Depending on the interval length W, theoptimal bandwidth w ∗ varies. We suggest selecting aninterval length that scales with the optimal bandwidthas γ −1 w∗ . The parameter γ regulates the interval length
for local optimization: With small γ ( 1), the fixedbandwidth is optimized within a long interval; Withlarge γ (∼ 1), the fixed bandwidth is optimized within ashort interval. The interval length and fixed bandwidth,γγ
selected at time t, are denoted as Wt and w̄t .
γ
The locally optimized bandwidth w̄t is repeatedly
obtained for different t(∈ [a, b ]). Because the intervals

J Comput Neurosci (2010) 29:171–182

175

overlap, we adopt the Nadaraya-Watson kernel regresγ
sion (Nadaraya 1964; Watson 1964) of w̄t as a local
bandwidth at time t:


γ
t−s γ
t−s
ρW
ρW
(17)
wt =
γ w̄s ds
γ ds.
s

s

γ

The variable bandwidth wt obtained from the samedata, but with different γ , exhibits different degrees
of smoothness: With small γ ( 1), the variable bandwidth fluctuates slightly; With large γ (∼ 1), the variablebandwidth fluctuates significantly. The parameter γ isthus a smoothing parameter for the variable bandwidth.
Similar to the fixed bandwidth, the goodness-of-fit of
the variable bandwidth can be estimated from the data.
The cost function for the variable bandwidth selectedwith γ is obtained as b2 Cˆn (γ ) =λ̂2t dt − 2
k γ ti − t j ,(18)n i= j wtia
where λ̂t = xt−s kwtγ (s) ds is an estimated rate, with theγvariable bandwidth wt . The integral is calculated numerically. With the stiffness constant γ ∗ that minimizes
Eq. (18), local optimization is performed in an ideal

3 Results
3.1 Comparison of the fixed and variable kernel
methods
By using spikes sampled from an inhomogeneous
Poisson point process, we examined the efficiency of
the kernel methods in estimating the underlying rate.
We also used a sequence obtained by superimposing
ten non-Poissonian (gamma) sequences (Shimokawa
and Shinomoto 2009), but there was practically nosignificant difference in the rate estimation from the
Poissonian sequence.
Figure 1 displays the result of the fixed kernel method based on the Gauss density function. The kernel
bandwidth selected by Algorithm 1 applies a reasonable filtering to the set of spike sequences. Figure 1(d)shows that a cost function, Eq. (10), estimated from the
spike data is similar to the original MISE, Eq. (4), which
was computed using the knowledge of the underlying

Spike Rate [spike/s](a)t

Underlying rate

20
0

(b)

1

2

Spike Sequences(d)Cn (w )

Cost function-100

Theoretical cost function

Empirical cost function

Too small w
0

(c)

1

2
-140

ˆt

Too small w

40
Too large w

20
Spike Rate [spike/s]

Fig. 1 Fixed kernel density
estimation. (a) The
underlying spike rate λt of the
Poisson point process. (b) 20
spike sequences sampledfrom the underlying rate, and(c) Kernel rate estimates
made with three types of
bandwidth: too small,
optimal, and too large. The
gray area indicates theunderlying spike rate. (d) The
cost function for bandwidth
w. Solid line is the estimatedcost function, Eq. (10),computed from the spike
data; The dashed line is theexact cost function, Eq. (5),directly computed by using
the known underlying rateinterval length. The method for optimizing the variable
kernel bandwidth is summarized in Algorithm 2.

-180

ˆt

Optimized w*

Optimized w *

20

0

0.1

0.2
Bandwidth w [s]ˆt
Too large w
20
0

1
Time [s]

2

0.3

0.4

176

J Comput Neurosci (2010) 29:171–182

rate. This demonstrates that MISE optimization can, in
practice, be carried out by our method, even without
knowing the underlying rate.
Figure 2(a) demonstrates how the rate estimationis altered by replacing the fixed kernel method withthe variable kernel method (Algorithm 2), for identical spike data (Fig. 1(b)). The Gauss weight function
is used to obtain a smooth variable bandwidth. The
manner in which the optimized bandwidth varies in thetime axis is shown in Fig. 2(b): the bandwidth is shortin a moment of sharp activation, and is long in the period of smooth rate modulation. Eventually, the sharp
activation is grasped minutely and slow modulation isexpressed without spurious fluctuations. The stiffness
constant γ for the bandwidth variation is selected byminimizing the cost function, as shown in Fig. 2(c).

3.2 Comparison with established density estimation
methods
We wish to examine the fitting performance of the
fixed and variable kernel methods in comparison withestablished density estimation methods, by paying attention to their aptitudes for either continuous or dis-continuous rate processes. Figure 3(a) shows the results
for sinusoidal and sawtooth rate processes, as samples of continuous and discontinuous processes, respectively. We also examined triangular and rectangularrate processes as different samples of continuous anddiscontinuous processes, but the results were similar.
The goodness-of-fit of the density estimate to the underlying rate is evaluated in terms of integrated squared
error (ISE) between them.
The established density estimation methods examined for comparison are the histogram (Shimazaki
and Shinomoto 2007), Abramson’s adaptive kernel(Abramson 1982), Locfit (Loader 1999b), and Bayesian
adaptive regression splines (BARS) (DiMatteo et al.
2001; Kass et al. 2003) methods, whose details are
summarized below.
A histogram method, which is often called a peristimulus time histogram (PSTH) in neurophysiological literature, is the most basic method for estimating
the spike rate. To optimize the histogram, we used amethod proposed for selecting the bin width based on
the MISE principle (Shimazaki and Shinomoto 2007).
Abramson’s adaptive kernel method (Abramson
1982)
uses the sample point kernel estimate λ̂t =k(t− ti ), in which the bandwidths are adaptedwtii

Spike Rate [spike/s](a)Underlying Rate

30

VARIABLE KERNELFIXED KERNEL

20

Cost function

Cn ( γ )Empirical cost funtion
Theoretical cost funtion

10

0

1

2

-185

Time [s](b)Bandwidth [s](c)Optimized Bandwidth
Optimized γ∗

VARIABLE (γ∗ = 0.8)
FIXED

0.2

-190

0.1

0
0

1

0.2

0.4
0.6
0.8
Stiffness constant γ

12Time [s]Fig. 2 Variable kernel density estimation. (a) Kernel rate estimates. The solid and dashed lines are rate estimates made by thevariable and fixed kernel methods for the spike data of Fig. 1(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 solidline is the cost function for the bandwidth stiffness constant γ ,Eq. (18), estimated from the spike data; the dashed line is thecost function computed from the known underlying rate

J Comput Neurosci (2010) 29:171–182

(a)

40

177

(b)

Sinusoid

400

40

Sawtooth2001Time [s]2

Underlying RateHISTOGRAM
FIXED KERNELVARIABLE KERNEL

Integrated Squared Error (Sawtooth)

Spike Rate [spikes/s]

20

HISTOGRAM
FIXED KERNELVARIABLE KERNELABRAMSON'SLOCFIT
BARS

300

200

100
ABRAMSON'SLOCFIT
BARS

0

50

100

150

200

250

Integrated Squared Error (Sinusoid)

Fig. 3 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 (grayarea), 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 intheir goodness-of-fit, based on the integrated squared error (ISE)
between the underlying and estimated rate. The abscissa andthe 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 plottedfor each method

at the sample points. Scaling the bandwidths as wti =w (g/λ̂ti )1/2 was suggested, where w is a pilot bandwidth,
g = ( i λ̂ti )1/N , and λ̂t is a fixed kernel estimate with w.
Abramson’s method is a two-stage method, in which
the pilot bandwidth needs to be selected beforehand.
Here, the pilot bandwidth is selected using the fixedkernel optimization method developed in this study.
The Locfit algorithm developed by Loader (1999b)fits a polynomial to a log-density function under the
principle of maximizing a locally defined likelihood.
We examined the automatic choice of the adaptive
bandwidth of the local likelihood, and found that thedefault fixed method yielded a significantly better fit.
We used a nearest neighbor based bandwidth method,with a parameter covering 20% of the data.
The BARS (DiMatteo et al. 2001; Kass et al. 2003)
is a spline-based adaptive regression method on an
exponential family response model, including a Poissoncount distribution. The rate estimated with the BARS
is the expected splines computed from the posterior
distribution on the knot number and locations witha Markov chain Monte Carlo method. The BARS is,thus, capable of smoothing a noisy histogram without
missing abrupt changes. To create an initial histogram,
we used 4 [ms] bin width, which is small enough to
examine rapid changes in the firing rate.

Figure 3(a) displays the density profiles of the six different methods estimated from an identical set of spiketrains (n = 10) that are numerically sampled from a sinusoidal or sawtooth underlying rate (2 [s]). Figure 3(b)summarizes the goodness-of-fit of the six methods tothe sinusoidal and sawtooth rates (10 [s]) by averaging
over 20 realizations of samples.
For the sinusoidal rate function, representing continuously varying rate processes, the BARS is mostefficient in terms of ISE performance. For the sawtooth rate function, representing discontinuous nonstationary rate processes, the variable kernel estimation
developed here is the most efficient in grasping abrupt
rate changes. The histogram method is always inferior
to the other five methods in terms of ISE performance,
due to the jagged nature of the piecewise constantfunction.
3.3 Application to experimental data
We examine, here, the fixed and variable kernel methods in their applicability to real biological data. In
particular, the kernel methods are applied to the spikedata of an MT neuron responding to a random dot
stimulus (Britten et al. 2004). The rates estimated fromn = 1, 10, and 30 experimental trials are shown in178

J Comput Neurosci (2010) 29:171–182

Fig. 4. Fine details of rate modulation are revealed
as we increase the sample size (Bair and Koch 1996).
The fixed kernel method tends to choose narrower
bandwidths, while the variable kernel method tends tochoose wider bandwidths in the periods in which spikes
are not abundant.
The performance of the rate estimation methods
is cross-validated. The bandwidth, denoted as wt for(a)both fixed and variable, is obtained with a training
data set of n trials. The error is evaluated by computing the cost function, Eq. (18), in a cross-validatory
manner:Ĉn (wt ) = ba

λ̂2t dt −



2 
t
,
k
−
t
w
i
j
n2 i= j ti

(19)

n=1

100

FIXED KERNEL50

VARIABLE KERNEL

0.4
0.2
0

Optimized Bandwidth

(d)

(b)

n=10

The difference in cross-validated MISE

100

100
50

50
0.04
0.02
0

0

-50

(c)

n=30

100
-100

50

1

5

10

20

30

The number of spike sequences, n

0.02
0.01
0

0

1
Time [s]2Fig. 4 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. 2004)). (a–c): Analyses of n = 1, 10,
and 30 spike trains; (top) Spike rates [spikes/s] estimated withthe fixed and variable kernel methods are represented by thegray area and solid line, respectively; (middle) optimized fixedand variable bandwidths [s] are represented by dashed and solidlines, respectively; (bottom) A raster plot of the spike data usedin the estimation. (d) Comparison of the two kernel methods.
Bars represent the difference between the cross-validated costfunctions, Eq. (19), of the fixed and variable kernel methods (fixless 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(ntotal = 60) are divided into ntotal /n blocks, each composedof n(= 1, 5, 20, 30) spike sequences. A bandwidth was selectedusing spike data of a training block. The cross-validated costfunctions, Eq. (19), for the selected bandwidth are computed
using the ntotal /n − 1 leftover test blocks, and their average iscomputed. The cost function is repeatedly obtained ntotal /n-times
by changing the training block. The mean and standard deviation,
computed from ntotal /n samples, are displayed

J Comput Neurosci (2010) 29:171–182

where the test spike times {ti } are obtained from nspike sequences in the leftovers, and λ̂t = n1 i kwt (t−ti ).
Figure 4(d) shows the performance improvements by
the variable bandwidth over the fixed bandwidth, as
evaluated by Eq. (19). The fixed and variable kernel
methods perform better for smaller and larger sizesof data, respectively. In addition, we compared thefixed kernel method and the BARS by cross-validating
the log-likelihood of a Poisson process with the rateestimated using the two methods. The difference inthe log-likelihoods was not statistically significant forsmall samples (n = 1, 5 and 10), while the fixed kernel
method fitted better to the spike data with larger samples (n = 20 and 30).

4 Discussion
In this study, we developed methods for selecting thekernel bandwidth in the spike rate estimation basedon the MISE minimization principle. In addition to theprinciple of optimizing a fixed bandwidth, we further
considered selecting the bandwidth locally in time, assuming a non-stationary rate modulation.
We tested the efficiency of our methods using
spike sequences numerically sampled from a given rate
(Figs. 1 and 2). Various density estimators constructed
on different optimization principles were compared intheir goodness-of-fit to the underlying rate (Fig. 3).
There is in fact no oracle that selects one among various
optimization principles, such as MISE minimization or
likelihood maximization. Practically, reasonable principles render similar detectability for rate modulation;
the kernel methods based on MISE were roughly comparable to the Locfit based on likelihood maximizationin their performances. The difference of the performances is not due to the choice of principles, but rather
due to techniques; kernel and histogram methods leadto completely different results under the same MISE
minimization principle (Fig. 3(b)). Among the smoothrate estimators, the BARS was good at representing
continuously varying rate, while the variable kernel
method was good at grasping abrupt changes in the rate
process (Fig. 3(b)).
We also examined the performance of our methods
in application to neuronal spike sequences by crossvalidating with the data (Fig. 4). The result demonstrated that the fixed kernel method performed well insmall samples. We refer to Cunningham et al. (2008)
for a result on the superior fitting performance of a179

fixed kernel to small samples in comparison with the
Locfit and BARS, as well as the Gaussian processsmoother (Cunningham et al. 2008; Smith and Brown
2003; Koyama and Shinomoto 2005). The adaptive
methods, however, have the potential to outperformthe fixed method with larger samples derived from a
non-stationary rate profile (See also Endres et al. 2008
for comparisons of their adaptive histogram with thefixed histogram and kernel method). The result in Fig. 4
confirmed the utility of our variable kernel method forlarger samples of neuronal spikes.
We derived the optimization methods under the
Poisson assumption, so that spikes are randomly drawn
from a given rate. If one wishes to estimate spike rateof a single or a few sequences that contain strongly correlated spikes, it is desirable to utilize the information
as to non-Poisson nature of a spike train (Cunningham
et al. 2008). Note that a non-Poisson spike train may bedually interpreted, as being derived either irregularly
from a constant rate, or regularly from a fluctuating
rate (Koyama and Shinomoto 2005; Shinomoto andKoyama 2007). However, a sequence obtained by superimposing many spike trains is approximated as a
Poisson process (Cox 1962; Snyder 1975; Daley andVere-Jones 1988; Kass et al. 2005), for which dual
interpretation does not occur. Thus the kernel methods
developed in this paper are valid for the superimposedsequence, and serve as the peristimulus density estimator for spike trains aligned at the onset or offset of thestimulus.
Kernel smoother is a classical method for estimating
the firing rate, as popular as the histogram method. We
have shown in this paper that the classical kernel methods perform well in the goodness-of-fit to the underlying rate. They are not only superior to the histogram
method, but also comparable to modern sophisticatedmethods, such as the Locfit and BARS. In particular,
the variable kernel method outperformed competing
methods in representing abrupt changes in the spike
rate, which we often encounter in neuroscience. Given
simplicity and familiarity, the kernel smoother can stillbe the most useful in analyzing the spike data, provided that the bandwidth is chosen appropriately as
instructed in this paper.
Acknowledgements We thank M. Nawrot, S. Koyama, D.
Endres for valuable discussions, and the Diesmann Unit forproviding the computing environment. We also acknowledge
K. H. Britten, M. N. Shadlen, W. T. Newsome, and J. A.
Movshon, who made their data available to the public, and
W. Bair for hosting the Neural Signal Archive. This study is
supported in part by a Research Fellowship of the Japan Society

180

J Comput Neurosci (2010) 29:171–182

for the Promotion of Science for Young Scientists to HS and
Grants-in-Aid for Scientific Research to SS from the MEXT
Japan (20300083, 20020012).
Open Access This article is distributed under the terms of theCreative Commons Attribution Noncommercial License which
permits any noncommercial use, distribution, and reproduction
in any medium, provided the original author(s) and source arecredited.

Appendix: Cost functions of the Gauss kernel function
In this appendix, we derive definite MISE optimization
algorithms we developed in the body of the paper withthe particular Gauss density function, Eq. (3).

z 2where erf (z) = √2π 0 e−t dt. A simplified equation is
obtained by evaluating the MISE in an unboundeddomain: a → −∞ and b → +∞. Using erf (±∞) = ±1,
we obtain(ti −t j )2
1
ψw ti , t j ≈ √e− 4w2 .
π 2w

(24)

Using Eq. (24) in Eq. (22), we obtain a formula for
selecting the bandwidth of the Gauss kernel function,√
2 π n2 Ĉn (w):2√ − (ti −t j )2N
2  − (ti −t2j )2+
e 4w − 2 2e 2w.
w
w i< j(25)

A.1 A cost function for a fixed bandwidth
A.2 A local cost function for a variable bandwidth
The estimated cost function is obtained, as in Eq. (10):n2 Ĉn (w) =ψw ti , t j − 2
kw ti − t j ,
i= j

i, j

The local cost function is obtained, as in Eq. (15):n2 Ĉnt (w, W) =where, from Eq. (11), bkw (t−ti ) kw t−t j dt.
ψw ti , t j =

i= j

u−t
kw (u − ti )kw (u − t j) ρW
du.

For the summations in the local cost function, Eq. (15),
we have equalities:
ti −t
kw ti − t j ρW
=

i= j

i, j



ψw (ti , ti ) + 2iψw ti , t j .

+2





kw ti − t j



t −t
ti −t
ρW
+ ρWj
,

i< j

(26)

(21)

i< j

By plugging Eqs. (20) and (21) into Eq. (10), the costfunction is simplified as .
n2 Ĉn (w) =ψw (ti , ti )i

ti −t
kw ti − t j ρW,tti , t j =
ψw,W

i< j

ψw ti , t j =where, from Eq. (16),

A symmetric kernel function, including the Gauss function, is invariant to exchange of ti and t j when computing kw (ti − t j). In addition, the correlation of the kernel
function, Eq. (11), is symmetric with respect to ti and t j.
Hence, we obtain the following relationships

kw ti − t j = 2
kw ti − t j ,
(20)


t
ti , t j − 2
ψw,W

i, j

a

i= j



ψw ti , t j − 2kw ti − t j



.

(22)

i< j

For the Gauss kernel function, Eq. (3), with bandwidth w, Eq. (11) becomes(ti −t j )2
1
ψw ti , t j = √e− 4w2π4w
2b − ti − t j
2a − ti − t j
× erf
− erf
,
2w
2w

(23)



t
ψw,W
ti , t j =



i, j

t
ψw,W
(ti , ti ) + 2

i



t
ψw,W
ti , t j .

i< j

(27)
The first equation holds for a symmetric kernel. The
second equation is derived because Eq. (16) is invariant
to an exchange of ti and t j. Using Eqs. (26) and (27),
Eq. (15) can be computed asn2 Ĉnt (w, W) =
+2





t
ψw,W
ti , t jitti , t j − kw ti − t j
ψw,W


!
t −t
ti −t
ρW
+ ρWj
.

i< j

(28)

J Comput Neurosci (2010) 29:171–182

181

For the Gauss kernel function and the Gauss weightfunction with bandwidth w and W respectively, Eq. (16)is calculated astψw,Wti , t j =


×

"

1
1
√
2π w2 2π W

(u − ti )2 + u − t j
exp −
2w 2

2

#
(u − t)2−
du. (29)
2W 2

By completing the square with respect to u, the exponent in the above equation is written as2
w2 +2W 2
τ W 2 + (ti −t) w2u+
2w2 W 2
w2 +2W 2(t − ti )2 + (t − t j)2 w2 + (ti − t j)2 W 2−.
2w 2 w 2 +2W 2
Using the formula
tained ase−Au du =
2

$π
A

(30)

, Eq. (16) is ob-

1√
2π w w 2 +2W 2
" 
#
 2
(t−ti )2 +(t−t j)2 w +(ti −t j)2 W 2
× exp −. (31)
2w 2 w2 +2W 2

ψ tw,W ti , t j =References
Abeles, M. (1982). Quantification, smoothing, and confidencelimits for single-units histograms. Journal of Neuroscience
Methods, 5(4), 317–325.
Abramson, I. (1982). On bandwidth variation in kernel estimatesa square root law. The Annals of Statistics, 10(4), 1217–1223.
Adrian, E. (1928). The basis of sensation: The action of the senseorgans. New York: W.W. Norton.
Bair, W., & Koch, C. (1996). Temporal precision of spike trains in
extrastriate cortex of the behaving macaque monkey. NeuralComputation, 8(6), 1185–1202.
Bowman, A. W. (1984). An alternative method of crossvalidation for the smoothing of density estimates. Biometrika, 71(2), 353.
Breiman, L., Meisel, W., & Purcell, E. (1977). Variable kernel estimates of multivariate densities. Technometrics, 19,
135–144.
Brewer, M. J. (2004). A Bayesian model for local smoothingin kernel density estimation. Statistics and Computing, 10,
299–309.
Britten, K. H., Shadlen, M. N., Newsome, W. T., & Movshon,
J. A. (2004). Responses of single neurons in macaque
mt/v5 as a function of motion coherence in stochastic
dot stimuli. The Neural Signal Archive. nsa2004.1. http://
www.neuralsignal.org.

Cox, R. D. (1962). Renewal theory. London: Wiley.
Cunningham, J., Yu, B., Shenoy, K., Sahani, M., Platt, J., Koller,D., et al. (2008). Inferring neural firing rates from spike trainsusing Gaussian processes. Advances in Neural Information
Processing Systems, 20, 329–336.
Daley, D., & Vere-Jones, D. (1988). An introduction to the theoryof point processes. New York: Springer.
Devroye, L., & Lugosi, G. (2000). Variable kernel estimates: On
the impossibility of tuning the parameters. In E. Giné, D.
Mason, & J. A. Wellner (Eds.), High dimensional probability
II (pp. 405–442). Boston: Birkhauser.
DiMatteo, I., Genovese, C. R., & Kass, R. E. (2001). Bayesian
curve-fitting with free-knot splines. Biometrika, 88(4),1055–1071.
Endres, D., Oram, M., Schindelin, J., & Foldiak, P. (2008).
Bayesian binning beats approximate alternatives: Estimating
peristimulus time histograms. Advances in Neural Information Processing Systems, 20, 393–400.
Fan, J., Hall, P., Martin, M. A., & Patil, P. (1996). On local
smoothing of nonparametric curve estimators. Journal of theAmerican Statistical Association, 91, 258–266.
Gerstein, G. L., & Kiang, N. Y. S. (1960). An approach to thequantitative analysis of electrophysiological data from single
neurons. Biophysical Journal, 1(1), 15–28.
Hall, P., & Schucany, W. R. (1989). A local cross-validation algorithm. Statistics & Probability Letters, 8(2), 109–117.
Jones, M., Marron, J., & Sheather, S. (1996). A brief survey of
bandwidth selection for density estimation. Journal of theAmerican Statistical Association, 91(433), 401–407.
Kass, R. E., Ventura, V., & Brown, E. N. (2005). Statistical issues
in the analysis of neuronal data. Journal of Neurophysiology,
94(1), 8–25.
Kass, R. E., Ventura, V., & Cai, C. (2003). Statistical smoothingof neuronal data. Network-Computation in Neural Systems,
14(1), 5–15.
Koyama, S., & Shinomoto, S. (2005). Empirical Bayes interpretations of random point events. Journal of Physics AMathematical and General, 38, 531–537.
Loader, C. (1999a). Bandwidth selection: Classical or plug-in?
The Annals of Statistics, 27(2), 415–438.
Loader, C. (1999b). Local regression and likelihood. New York:
Springer.
Loftsgaarden, D. O., & Quesenberry, C. P. (1965). A nonparametric estimate of a multivariate density function. The Annals of Mathematical Statistics, 36, 1049–1051.
Nadaraya, E. A. (1964). On estimating regression. Theory ofProbability and its Applications, 9(1), 141–142.
Nawrot, M., Aertsen, A., & Rotter, S. (1999). Single-trial estimation of neuronal firing rates: From single-neuron spiketrains to population activity. Journal of Neuroscience Methods, 94(1), 81–92.
Parzen, E. (1962). Estimation of a probability density-function
and mode. The Annals of Mathematical Statistics, 33(3),1065.
Richmond, B. J., Optican, L. M., & Spitzer, H. (1990). Temporal encoding of two-dimensional patterns by single units inprimate primary visual cortex. i. stimulus-response relations.
Journal of Neurophysiology, 64(2), 351–369.
Rosenblatt, M. (1956). Remarks on some nonparametric estimates of a density-function. The Annals of Mathematical
Statistics, 27(3), 832–837.
Rudemo, M. (1982). Empirical choice of histograms and kernel
density estimators. Scandinavian Journal of Statistics, 9(2),65–78.

182
Sain, S. R. (2002). Multivariate locally adaptive density estimation. Computational Statistics & Data Analysis, 39, 165–186.
Sain, S., & Scott, D. (1996). On locally adaptive density estimation. Journal of the American Statistical Association, 91(436),
1525–1534.
Sain, S., & Scott, D. (2002). Zero-bias locally adaptive density estimators. Scandinavian Journal of Statistics, 29(3),441–460.
Sanderson, A. (1980). Adaptive filtering of neuronal spike train
data. IEEE Transactions on Biomedical Engineering, 27,
271–274.
Scott, D. W. (1992). Multivariate density estimation: Theory, practice, and visualization. New York: Wiley-Interscience.
Scott, D. W., & Terrell, G. R. (1987). Biased and unbiased crossvalidation in density estimation. Journal of the American
Statistical Association, 82, 1131–1146.
Shimazaki, H., & Shinomoto, S. (2007). A method for selecting
the bin size of a time histogram. Neural Computation, 19(6),1503–1527.
Shimokawa, T., & Shinomoto, S. (2009). Estimating instantaneous irregularity of neuronal firing. Neural Computation,
21(7), 1931–1951.

J Comput Neurosci (2010) 29:171–182
Shinomoto, S., & Koyama, S. (2007). A solution to the controversy between rate and temporal coding. Statistics in Medicine, 26, 4032–4038.
Shinomoto, S., Kim, H., Shimokawa, T., Matsuno, N., Funahashi,S., Shima, K., et al. (2009). Relating neuronal firing patterns
to functional differentiation of cerebral cortex. PLoS Computational Biology, 5, e1000433.
Shinomoto, S., Miyazaki, Y., Tamura, H., & Fujita, I. (2005)
Regional and laminar differences in in vivo firing patterns ofprimate cortical neurons. Journal of Neurophysiology, 94(1),567–575.
Shinomoto, S., Shima, K., & Tanji, J. (2003). Differences in spiking patterns among cortical neurons. Neural Computation,
15(12), 2823–2842.
Silverman, B. W. (1986). Density estimation for statistics and data
analysis. London: Chapman & Hall.
Smith, A. C., & Brown, E. N. (2003). Estimating a state-spacemodel from point process observations. Neural Computation, 15(5), 965–991.
Snyder, D. (1975). Random point processes. New York: Wiley.
Watson, G. S. (1964). Smooth regression analysis. Sankhya: The
Indian Journal of Statistics, Series A, 26(4), 359–372.