Neural Coding 2007 Abstract p120-123, Montevideo, UruguayKernel Width Optimization in the Spike-rate Estimation
Hideaki Shimazaki
Department of Physics, Kyoto UniversityKyoto 606-8502, Japan
shimazaki@ton.scphys.kyoto-u.ac.jp
Shigeru Shinomoto
Department of Physics, Kyoto University
Kyoto 606-8502, Japanshinomoto@scphys.kyoto-u.ac.jp

ABSTRACTA classical tool for estimating the neuronal spike rate is a peri-stimulus time histogram (PSTH)
constructed from spike sequences aligned at the onset of a stimulus repeatedly applied to ananimal. We have recently established a method for selecting the bin size, so that the PSTH
best represents the unknown underlying rate [1]. The goodness of the fit we adopted as theoptimization principle is minimizing the mean integrated squared error (MISE) between theunderlying rate λt and the PSTH λ̂t ,Z b
MISE =aE(λt − λ̂t )2 dt,(1)where E refers to the expectation with respect to the spike generation process under a given
time-dependent rate λt . The method allows us to minimize the MISE from spike count statistics
alone, without knowing the underlying rate.
In this contribution, we consider a kernel rate estimator as λ̂t and suggest a method to select thewidth of a kernel under the MISE criterion. Generally, the cross-validation method is applicableto the least squares minimization [2, 3]. Here, we estimate the MISE fully utilizing the Poissonian
nature of spikes, as we have done in the PSTH optimization. The Poissonian assumption holds
in the limit of large number of trials, because spikes repeatedly recorded from a single neuron
under identical experimental conditions are in the majority mutually independent.
For a small number of spike sequences generated from modestly fluctuating underlying rate,
the optimal kernel width may become comparable to the observation period. This phenomenon,
similar to what we have observed in the PSTH optimization [4], indicates that more experimental
trials are needed if one wishes to uncover the time-dependent rate. We also construct a method
for estimating the number of additional experimental trials needed to analyze the data with theresolution we deem sufficient.
Selection of the Kernel Width We consider independently and identically obtained n spikesequences, which contain N spikes as a whole. Due to the general limit theorem for the sumof independent point processes, a superposition of the spike sequences can be regarded as being
drawn from
P a time-dependent Poisson point process. We define the superimposed sequence as
xt = n−1 N
i=1 δ (t − ti ), where ti is the timing of the ith spike. δ(t) is the Dirac delta function.
An estimator
λ̂t is constructed by smoothing the point process by a kernel kw (t) of the bandwidthRw: λ̂t = kw (t − s)xs ds.
We wish to obtain a kernel function that minimizes the MISE (Eq. 1). The integrand of theMISE is decomposed into three parts: λ2t − 2λt E λ̂t + E λ̂2t . Since the first component does notdepend on the choice of a kernel, we subtract it from the MISE and define a cost function as a

Table 1: A Method for Selecting a Bandwidth of a Kernel Function
(i)
(ii)

(iii)

Superimpose all the n spike sequences. Obtain a series of spike times {ti }N
i=1in [a, b]. N is the total number of spikes.
Compute the cost function of a kernel kw (t) as
4 X
1 X
Ĉn (w) = − 2kw (ti − tj ) + 2
ψw,a,b (ti − tj ) ,nni<ji,jRb
where ψw,a,b (t) ≡ a kw (s) kw (s + t) ds is the correlation function (∗) .
Repeat ii while changing w to search for w∗ that minimizes Ĉn (w).
³´t2− 2w
2 ,´© ¡³¢
¡ 2a+t ¢ª1t2
ψw,a,b (t) = √π4w
−erf.
exp − 4w
erf 2b+t22w
2w(∗) For a Gaussian kernel k1w (t) = √2πw exp

function of the kernel, or its bandwidth w,Z bCn (w) ≡ MISE −a

λ2t dt = −2Z baZ bλt E λ̂t dt +a

E λ̂2t dt.

(2)

From a general decomposition rule of a covariance of two random variables, we obtain therelation,Z baZ b
λt E λ̂t dt =Z bE(xt − Ext )(λ̂t − E λ̂t ) dt
Z b
Z bkw (0)Ext dt.
=Ext λ̂t dt −naaaE[xt λ̂t ] dt −a(3)To obtain the last equality, we used the assumption that the spikes are independent each other(an assumption of a Poisson point process). Hence from sample sequences, the cost function isestimated asZZ b
Z b
2kw (0) b
Ĉn (w) =xt dt − 2xt λ̂t dt +λ̂2t dt
naaaN NN N
2kw (0)
2 XX
1 XX=N−k(t−t)+ψw,a,b (ti − tj ) ,(4)w i
jn2n2n2i=1 j=1i=1 j=1

where ψw,a,b (t) is given byZ b
ψw,a,b (t) ≡akw (s) kw (s + t) ds.

(5)

We summarize the method for selecting the width of a symmetric kernel function as Table 1.
Figure 1 displays application of our method to the spike data of an MT neuron responding toa random dot stimulus [5]. In this example, we changed the number of spike sequences used toestimate the time-dependent rate.
Extrapolation of the cost function Here we provide the method to estimate how many
additional experimental trials are needed to make an optimized rate estimate with the required
resolution. Suppose we have n spike sequences. The superposition of n spike sequences is(n)(n)denoted as xt . The kernel estimator is denoted as λ̂t .

B80

w*=90 [ms]n=1w*=13 [ms]n=5w*=6 [ms]n=30400ACost Function

C n (w)

C

80
40

Optimal Bandwidth
w*=90 [ms] (n=1)-200

0-400Dw*=13 [ms] (n=5)

80
40-600w*=6 [ms] (n=30)
0

Bandwidth w [s]





Time [s]Figure 1: Optimization of a kernel¡ 2estimate¢ √for spike-rate of an MT neuron (j024 in nsa2004.12[5]) with a Gaussian kernel (exp −t /2w / 2πw), according to the method in Table 1. A: Cost
functions of n = 1, 5, and 30 spike sequences. B, C, and D: (Top) Optimized time-dependent
rates. (Bottom) Spike sequences.
We wish to compute the cost function for m spike sequences. It is easily proved that the costfunction for m spike sequences is written asµ¶Z b Z
11Cm (w) =−dt kw (t − s)2 Ex(n)(6)
s ds + Cn (w),m naR(m)(m)
where Cn (w) is given by Eq. 2. To obtain Eq. 6, we used a relation, E(λ̂t − E λ̂t )2 dt =RR(m)(m)(n)
m−1 dt kw (t − s)2 Exs ds, and replaced Exs with Exs . Hence the cost function for m
spike sequences is estimated from a sample of n spike sequences asµĈm [w | n] =11−m n¶N1Xni=1Z b−tia−tikw (t)2 dt + Ĉn (w).

(7)

The extrapolation method is summarized as Table 2. Figure 2 demonstrates the application ofthe extrapolation method to the spike data [5]. In Figure 2-b, the optimal widths were estimatedwith the extrapolated cost functions computed from the first two spike sequences (Solid line).
They are close to the optimal kernel widths obtained by the method in Table 1 (Dots). With thisplot, experimentalists can estimate how many additional experimental trials must be performed
in order to achieve the resolution they deem sufficient.
Keywords: Optimization, Kernel Bandwidth, Rate Estimation, Spike Sequences

Table 2: A Method for Extrapolating the Cost Function(A)(B)

Compute the cost function Ĉn (w) as in Table 1.
Construct the extrapolated cost function:µ¶N
11 1X
Ĉm (w | n) =−ψw,a−ti ,b−ti (0) + Ĉn (w) .
m n n(C)

Repeat B while changing w to search for w∗ that minimizes Ĉm (w | n).

i=1

A Cn(w)

B w*
Cost function (n=2)
Extrapolation (m=30)0Extrapolation (m=60)

10-2-400

10-80000.02

0.04
Bandwidth w [s]Optimal Bandwidth
Extrapolation (n=2)

-3

100

11010
The number of sequences

2Figure 2: Application of the extrapolation method to spike data of an MT neuron (the data
as in Fig. 1 [5]). A: A cost function of the first two spike sequences (Solid). Extrapolated
cost functions for m = 30 and 60 computed from the first two spike sequences (Dashed andDotted). B: Estimated optimal kernel widths of n = 1 through 60 spike sequences, according to
the method in Table 1 (Dots). A solid line is the estimated optimal widths computed from thefirst 2 spike sequences, according to the method in Table 2.

References[1] Shimazaki, H. and Shinomoto, S. (2007) A method for selecting the bin size of a timehistogram, Neural Computation, 19: 1503-1527.
[2] Rudemo, M. (1982) Empirical choice of histograms and kernel density estimators, Scandinavian Journal of Statistics, 9: 65-78.
[3] Bowman, A. W. (1984) An alternative method of cross-validation for the smoothing of density
estimates, Biometrika, 71: 353-360.
[4] Koyama, S. and Shinomoto, S. (2004) Histogram bin width selection for time-dependent
Poisson processes, Journal of Physics A-Mathematical and General, 37: 7255-7265.
[5] Britten, K. H., Shadlen, M. N., Newsome, W. T., and Movshon, J. A. (2004). Responsesof single neurons in macaque MT/V5 as a function of motion coherence in stochastic dot
stimuli. Neural Signal Archive, nsa2004.1. http://www.neuralsignal.org.