LETTER

Communicated by Robert Kass

A Method for Selecting the Bin Size of a Time Histogram
Hideaki Shimazaki
shimazaki@ton.scphys.kyoto-u.ac.jp

Shigeru Shinomoto
shinomoto@scphys.kyoto-u.ac.jp
Department of Physics, Kyoto University, Kyoto 606-8502, JapanThe time histogram method is the most basic tool for capturing a timedependent rate of neuronal spikes. Generally in the neurophysiologicalliterature, the bin size that critically determines the goodness of the fitof the time histogram to the underlying spike rate has been subjectively
selected by individual researchers. Here, we propose a method for objectively selecting the bin size from the spike count statistics alone, sothat the resulting bar or line graph time histogram best represents theunknown underlying spike rate. For a small number of spike sequences
generated from a modestly fluctuating rate, the optimal bin size may
diverge, indicating that any time histogram is likely to capture a spurious rate. Given a paucity of data, the method presented here can nevertheless suggest how many experimental trials should be added in order to obtain a meaningful time-dependent histogram with the required
accuracy.
1 Introduction
Neurophysiological studies are based on the idea that information is transmitted between cortical neurons by spikes (Johnson, 1996; Dayan & Abbott,2001). A number of filtering algorithms have been proposed for estimating the instantaneous activity of an individual neuron or the joint activityof multiple neurons (DiMatteo, Genovese, & Kass, 2001; Wiener & Richmond, 2002; Sanger, 2002; Kass, Ventura, & Cai, 2003; Brockwell, Rojas, &Kass, 2004; Kass, Ventura, & Brown, 2005; Brown, Kass, & Mitra, 2004). The
most basic and frequently used tool for spike rate estimation is the time
histogram method. For instance, one aligns spike sequences at the onsetof stimuli repeatedly applied to an animal and describes the response of asingle neuron with a peristimulus time histogram (PSTH) or the responsesof multiple neurons with a joint PSTH (Adrian, 1928; Gerstein & Kiang,
1960; Gerstein & Perkel, 1969; Abeles, 1982).
The shape of a PSTH is largely dependent on the choice of the bin size.
With a bin size that is too large, one cannot represent the time-dependent
spike rate. And with a bin size that is too small, the time histogram
Neural Computation 19, 1503–1527 (2007)


C 2007 Massachusetts Institute of Technology

1504

H. Shimazaki and S. Shinomoto

fluctuates greatly, and one cannot discern the underlying spike rate. Thereis an appropriate bin size for each set of spike sequences, which is based on
the goodness of the fit of the PSTH to the underlying spike rate. For mostpreviously published PSTHs, however, the bin size has been subjectively
selected by the authors.
For data points distributed compactly, there are classical theories about
how the optimal bin size scales with the total number of data points n. It
was proven that the optimal bin size scales as n−1/3 with regard to the bargraph-density estimator (Révész, 1968; Scott, 1979). It was recently foundthat for two types of infinitely long spike sequences, whose rates fluctuate
either smoothly or jaggedly, the optimal bin sizes exhibit different scaling
relations with respect to the number of sequences, timescale, and amplitudeof rate modulation (Koyama & Shinomoto, 2004).
Though interesting, the scaling relations are valid only for a large amountof data and are of limited use in selecting a bin size. We devised a methodof selecting the bin size of a time histogram from the spike data. In thecourse of our study, we realized that a theory on the empirical choice ofthe histogram bin size for a probability density function was presented
by Rudemo (1982). Although applicable to a Poisson point process, thistheory appears to have rarely been used by neurophysiologists in analysesof PSTHs. In the actual procedure of neurophysiological experiments, thenumber of trials (spike sequences) plays an important role in determining
the resolution of a PSTH and thus in designing experiments. Therefore, itis preferable to have a theory that accords with the common protocol ofneurophysiological experiments in which a stimulus is repeated to extract
a signal from a neuron. Given a set of experimental data, we wish to not
only determine the optimal bin size, but also estimate how many moreexperimental trials should be performed in order to obtain a resolution wedeem sufficient.
For a small number of spike sequences derived from a modestly fluctuating rate, the estimated optimal bin size may diverge, implying that byconstructing a PSTH, it is likely that one obtains spurious results for the
spike rate estimation (Koyama & Shinomoto, 2004). Because a shortage ofdata underlies this divergence, one can carry out more experiments to obtaina reliable rate estimation. Our method can suggest how many sequences
should be added in order to obtain a meaningful time histogram with therequired accuracy. As an application of this method, we also show that thescaling relations of the optimal bin size that appears for a large number of
spike sequences can be examined from a relatively small amount of data.
The degree of the smoothness of an underlying rate process can be estimated by this method. In addition to a bar graph (piecewise constant) time
histogram, we also designed a method for creating a line graph (piecewise
linear) time histogram, which is superior to a bar graph in the goodness ofthe fit to the underlying spike rate and in comparing multiple responses todifferent stimulus conditions.

The Bin Size of a Time Histogram

1505

These empirical methods for the bin size selection for a bar and a line
graph histogram, estimation of the number of sequences required for thehistogram, and estimation of the scaling exponents of the optimal bin size
were corroborated by theoretical analysis derived for a generic stochastic
rate process. In the next section, we develop the bar graph (peristimulus)
time histogram (Bar-PSTH) method, which is the most frequently used
PSTH. In the appendix, we develop the line graph (peristimulus) time
histogram (Line-PSTH) method.
2 Optimization of the Bar Graph Time Histogram
We consider sequences of spikes repeatedly recorded from a single neuron
under identical experimental conditions. A recent analysis revealed thatin vivo spike trains are not simply random, but possess interspike interval distributions intrinsic and specific to individual neurons (Shinomoto,
Shima, & Tanji, 2003; Shinomoto, Miyazaki, Tamura, & Fujita, 2005). However, spikes accumulated from a large number of spike trains are in the
majority mutually independent and can be regarded as being derived froma time-dependent Poisson point process (Snyder, 1975; Daley & Vere-Jones,
1988; Kass et al., 2005).
It would be natural to assess the goodness of the fit of the estimator λ̂t
to the underlying spike rate λt over the total observation period T by themean integrated squared error (MISE),
1MISE ≡T T

E (λ̂t − λt )2 dt,

(2.1)

0

where E refers to the expectation over different realizations of point events,
given λt . We begin with a bar graph time histogram as λ̂t and explore amethod to select the bin size that minimizes the MISE. The difficulty of the
problem comes from the fact that the underlying spike rate λt is not known.
A bar graph time histogram is constructed simply by counting the number of spikes that belong to each bin of width . For an observation period T,
we obtain N = T/ intervals. The number of spikes accumulated fromall n sequences in the ith interval is counted as ki . The bar height at theith bin is given as ki /n. Figure 1 shows the schematic diagram for theconstruction of a bar graph time histogram.
Given a bin of width , the expected height of a bar graph for t ∈ [0, ]
is the time-averaged rate,

θ=

1


 
0

λt dt.

(2.2)

1506

H. Shimazaki and S. Shinomoto

A

0

B

^

C
0

Figure 1: Bar-PSTH. (A) An underlying spike rate, λt . The horizontal bars indicate the time-averaged rates θ for individual bins of width . (B) Sequences of
spikes derived from the underlying rate. (C) A time histogram for the sample
sequences of spikes. The estimated rate θ̂ is the total number of spikes k thatentered each bin, divided by the number of sequences n and the bin size .

The total number of spikes k from n spike sequences that enter this bin ofwidth  obeys the Poisson distribution,p(k | nθ ) =(nθ )k −nθ
e
,
k!

(2.3)

whose expectation is nθ . The unbiased estimator for θ is denoted as θ̂ =k/(n), which is the empirical height of the bar graph for t ∈ [0, ].
2.1 Selection of the Bin Size. By segmenting the total observation period T into N intervals of size , the MISE defined in equation 2.1 can be
rewritten as
1
MISE =


 
0


1 
E ( θ̂i − λt+(i−1) )2 dt,
N
N

(2.4)

i=1

where θ̂i ≡ ki /(n). Hereafter, we denote the average over the segmented
rate λt+(i−1) as an average over an ensemble of (segmented) rate functions

The Bin Size of a Time Histogram

1507{λt } defined on an interval of t ∈ [0, ], as
MISE =

1


 




E ( θ̂ − λt )2 dt.

(2.5)

0

The expectation E refers to the average over the spike count, or θ̂ = k/(n),given a rate function λt , or its mean value θ .
The MISE can be decomposed into two parts:
MISE = E(θ̂ − θ )2  +

1


 




(λt − θ )2 dt.

(2.6)

0

The first and second terms are, respectively, the stochastic fluctuation of theestimator θ̂ around the expected mean rate θ and the averaged temporalfluctuation of λt around its mean θ over an interval of length .
The second term of equation 2.6 can be decomposed further into twoparts:1 

(λt − θ  + θ  − θ )2 dt =

0

1


 

(λt − θ )2  dt − (θ − θ )2 .

0

(2.7)
The first term in equation 2.7 represents a mean squared fluctuation of theunderlying rate λt from the mean rate θ  and is independent of the bin size
, because
1


 

(λt − θ )2  dt =

0

1
T

 T

(λt − θ )2 dt.

(2.8)

0

We define a cost function by subtracting this term from the original MISE,
Cn () ≡ MISE −

1
T

 T

(λt − θ )2 dt

0= E(θ̂ − θ )2  − (θ − θ )2 .

(2.9)

This cost function corresponds to the risk function in Rudemo (1982, equation 2.3), obtained by direct decomposition of the MISE. The second term inequation 2.9 represents the temporal fluctuation of the expected mean rate
θ for individual intervals of period . As the expected mean rate θ is not
an observable quantity, we have to replace the fluctuation of the expected
mean rate with that of the observable estimator θ̂ . Using the decomposition

1508

H. Shimazaki and S. Shinomotorule for an unbiased estimator (E θ̂ = θ ),E(θ̂ − E θ̂ )2  = E(θ̂ − θ )2  + (θ − θ )2 ,
the cost function is transformed into 
Cn () = 2 E(θ̂ − θ )2 − E(θ̂ − E θ̂ )2 .

(2.10)(2.11)

Due to the assumed Poisson nature of the point process, the number of
spikes k counted in each bin obeys a Poisson distribution; the variance of
k is equal to the mean. For the estimated rate defined as θ̂ = k/(n), thisvariance-mean relation corresponds toE(θ̂ − θ )2 =

1
E θ̂.
n

(2.12)

By incorporating equation 2.12 into equation 2.11, the cost function is given
as a function of the estimator θ̂ ,Cn () =
2   E θ̂ − E(θ̂ − E θ̂)2 .
n

(2.13)The optimal bin size is obtained by minimizing the cost function Cn (), as∗ ≡ arg min Cn ().

(2.14)



By replacing the expectation with the sample spike count, the cost function(see equation 2.13) is converted into this useful recipe:
Algorithm 1: A Method for Bin Size Selection for a Bar-PSTH
i. Divide the observation period T into N bins of width , and countthe number of spikes ki from all n sequences that enter the ith bin.
ii. Construct the mean and variance of the number of spikes {ki } as
k̄ ≡

1 
1 
ki , and v ≡(ki − k̄)2 .
N
N
N

N

i=1

i=1

iii. Compute the cost function:
Cn () =2k̄ − v.
(n)2

iv. Repeat i through iii while changing the bin size  to search for ∗that minimizes Cn ().

The Bin Size of a Time Histogram1509

2.2 Extrapolation of the Cost Function. With the method developed in
the preceding section, we can determine the optimal bin size for a given setof experimental data. In this section, we develop a method to estimate how
the optimal bin size decreases when more experimental trials are added tothe data set: given n sequences, the method provides the cost function form(= n) sequences.
Assume that we are in possession of n spike sequences. The fluctuationof the expected mean rate (θ − θ )2  in equation 2.10 is replaced with theempirical fluctuation of the time histogram θ̂ n using the decomposition rule
for the unbiased estimator θ̂n satisfying E θ̂n = θ ,


 
 

E(θ̂n − E θ̂n )2 = E(θ̂n − θ )2 + (θ − θ )2 .

(2.15)The expected cost function for m sequences can thus be obtained by substituting the above equation into equation 2.9, yielding

 
 

Cm (|n) = E(θ̂m − θ )2 + E(θ̂n − θ )2 − E(θ̂n − E θ̂n )2 .

(2.16)

Using the variance-mean relation for a Poisson distribution, equation 2.12,
andE(θ̂m − θ )2 =1
1
E θ̂m =E θ̂n ,mm(2.17)

we obtainCm (|n) =

1
1
−
m n




1 
E θ̂n + Cn () ,(2.18)

where Cn () is the original cost function (see equation 2.13) computed
using the estimators θ̂n . By replacing the expectation with sample spikecount averages, the cost function for m sequences can be extrapolated withthis formula, using the sample mean k̄ and variance v of the numbers of
spikes, given n sequences and the bin size . The extrapolation method issummarized as algorithm 2:
Algorithm 2: A Method for Extrapolating the Cost Function for a
Bar-PSTH
A. Construct the extrapolated cost function,Cm (|n) =

1
1
−
m n



k̄
+ Cn (),
n2

1510

H. Shimazaki and S. Shinomotousing the sample mean k̄ and variance v of the number of spikes
obtained from n sequences of spikes. Cn () is the cost function
computed for n sequences of spikes with algorithm 1.
B. Search for ∗m that minimizes Cm (|n).
C. Repeat A and B while changing m, and plot 1/∗m versus 1/m
to search for the critical value 1/m = 1/n̂c above which 1/∗mpractically vanishes.
It may come to pass that the original cost function Cn () computed forn spike sequences does not have a minimum or has a minimum at a bin
size comparable to the observation period T. Because of the paucity of
data, one may consider carrying out more experiments to obtain a reliable
rate estimation. The critical number of sequences nc above which the costfunction has a finite bin size ∗ may be estimated in the following manner.
With a large , the cost function can be expanded asCn () ∼ µ

1
1
−
n nc



1
1
+u 2,



(2.19)

where we have introduced nc and u, which are independent of n. The
optimal bin size undergoes a phase transition from the vanishing 1/∗ forn < nc to a finite 1/∗ for n > nc . The inverse optimal bin size is expandedin the vicinity of nc as 1/∗ ∝ (1/n − 1/nc ). We can estimate the criticalˆ ∗mvalue n̂c (see Figure 3) by applying this asymptotic relation to the set of 
estimated from Cm (|n) for various values of m:


1
1
1
−
.
(2.20)∝
∗m
m n̂c
The minimum number of sequences required for the construction of a BarPSTH is estimated from n̂c . It should be noted that there are cases that the
optimal bin size exhibits a discontinuous divergence from a finite value.
Even in such cases, the plot of {1/m, 1/∗ } could be useful in exploring thediscontinuous transition from nonvanishing values of 1/∗ to practicallyvanishing values.
In the opposite extreme, with a sufficiently large number of spike sequences, our method selects a small bin size. It is known that the optimal
bin size exhibits a power law scaling with respect to the number of sequences n. The exponent of the scaling relation depends on the smoothnessof the underlying rate (Koyama & Shinomoto, 2004). Given a large numberof spike sequences, the method for extrapolating the cost function can also
be used to estimate the optimal bin size for a larger number of spike sequence m > n, and further estimate the scaling exponent representing thesmoothness of the underlying rate.

The Bin Size of a Time Histogram

1511

2.3 Theoretical Analysis on the Optimal Bin Size. To verify the empirical methods, we obtain the theoretical cost function of a Bar-PSTH directlyfrom a process with a known underlying rate. Note that this theoreticalcost function is not available in real experimental conditions in which theunderlying rate is not known.
The estimator θ̂ ≡ k/(n) is a uniformly minimum variance unbiased
estimator (UMVUE) of θ , which achieves the lower bound of the CramérRao inequality (Blahut, 1987; Cover & Thomas, 1991),

E(θ̂ − θ )2 = −

∞

k=0

p (k | nθ )∂ 2 log p (k | nθ )
∂θ 2

−1

=

θ
.
n

(2.21)

Inserting this into equation 2.9, the cost function is represented asθ − (θ − θ )2 n  
µ1=−
φ (t1 − t2 ) dt1 dt2 ,n 2 0 0

Cn () =(2.22)

where µ = θ  is the mean rate, and φ(t1 − t2 ) = (λt1 − µ)(λt2 − µ) is the
autocorrelation function of the rate fluctuation, λt − µ. To obtain the lastequation, we used


2
1 
(θ − θ )  =(λt − µ) dt 0  
1
(λt1 − µ)(λt2 − µ) dt1 dt2 .
= 2
 0 0
2

(2.23)

The cost function with a large bin size can be rewritten as
 
µ
1
− 2
( − |t|)φ(t) dtn  −
 ∞

1 ∞
1
µ
−
φ(t) dt + 2|t|φ(t) dt,∼n  −∞ −∞Cn () =(2.24)

based on the symmetry φ(t) = φ(−t) for a stationary process. Equation 2.24
can be identified with equation 2.19 with parameters
µ
φ(t)
dt
−∞

nc =  ∞(2.25)

1512

H. Shimazaki and S. Shinomoto

u=

 ∞
−∞

|t|φ(t) dt.

(2.26)For the process giving the correlation of the form φ (t) = σ 2 e −|t|/τ and the22
process giving a gaussian correlation φ(t) = σ 2 e −t /τ used in the simulationin the next √section, the critical numbers are obtained as nc = µ/2σ 2 τ and2
nc = µ/σ τ π , respectively.
Based on the same theoretical cost function, equation 2.22, we can also
derive scaling relations of the optimal bin size, which is achievable witha large number of spike sequences. With a small , the correlation of rate
fluctuation φ(t) can be expanded as φ (t) = φ (0) + φ  (0+ ) |t| + 12 φ  (0) t 2 + 
O |t|3 , and we obtain an expansion of equation 2.22 with respect to ,Cn () =1µ1− φ (0) − φ  (0+ )  − φ  (0) 2 + O(3 ).
n
3
12(2.27)The optimal bin size ∗ is obtained from dCn () /d = 0. For a rate thatfluctuates smoothly in time, the correlation function is a smooth function oft, resulting in φ  (0) = 0 due to the symmetry φ(t) = φ(−t). In this case, we
obtain the scaling relation6µ ∼ − φ (0) n∗

1/3.

(2.28)For a rate that fluctuates in a zigzag pattern, in which the correlationof rate fluctuation has a cusp at t = 0 (φ  (0+ ) < 0), we obtain the scaling
relation∗ ∼ −3µ
φ (0+ )n

1/2,

(2.29)

by ignoring the second-order term. Examples that exhibit this type of scaling are the Ornstein-Uhlenbeck process and the random telegraph process (Kubo, Toda, & Hashitsume, 1985; Gardiner, 1985).
The scaling relations, equation 2.28 and equation 2.29, are the generalization of the ones found in Koyama and Shinomoto (2004) for two specifictime-dependent Poisson processes.
3 Application of the Method to Spike Data
In this section, we apply the methods for the bar graph (peristimulus) time
histogram (Bar-PSTH) and the line graph (peristimulus) time histogram(Line-PSTH) to (1) spike sequences (point events) generated by the simulations of time-dependent Poisson processes and (2) spike sequences recorded

The Bin Size of a Time Histogram1513

from a neuron in area MT. The methods for a Line-PSTH are summarized
in the appendix.
3.1 Application to Simulation Data. We applied the method for estimating the optimal bin size of a Bar-PSTH to a set of spike sequences
derived from a time-dependent Poisson process. Figure 2A shows the empirical cost function Cn () computed from a set of n spike sequences. The
empirical cost function is in good agreement with the theoretical cost function computed from the mean spike rate µ and the correlation φ(t) of therate fluctuation, according to equation 2.22. In Figure 2D, a time histogramwith the optimal bin size is compared with those with nonoptimal bin sizes,
demonstrating the effectiveness of optimizing the bin size.
Algorithm 2 provides an estimate of how many sequences should be
added for the construction of a Bar-PSTH with the resolution we deem
sufficient. In this method, the cost function of m sequences is obtainedby modifying the original cost function, Cn (), computed from the spikecount statistics of n sequences of spikes. In subsequent applications of thisextrapolation method, the original cost function, Cn (), was obtained byaveraging over the initial partitioning positions.
We applied this extrapolation method for a Bar-PSTH to a set of spikesequences derived from the smoothly regulated Poisson process. Figure 3A
depicts the extrapolated cost function Cm (|n) for several values of m, computed from a given set of n(= 30) spike sequences. Figure 3B represents thedependence of an inverse optimal bin size 1/∗m on 1/m. The inverse optimal bin size, 1/∗m , stays near 0 for 1/m > 1/n̂c and departs from 0 linearlywith m for 1/m ≤ 1/ n̂c . By fitting the linear function, we estimated the critical number of sequences n̂c . In Figure 3C, the critical number of sequencesn̂c estimated from a smaller or larger n is compared to the theoretical valueof nc computed from equation 2.25. It is noteworthy that the estimated n̂capproximates the theoretical nc well even from a fairly small number of
spike sequences, with which the estimated optimal bin size diverges.
We also constructed a method for selecting the bin size of a Line-PSTH,
which is summarized as algorithm 3 in the appendix. Figure 4 compares the
optimal Bar-PSTH and the optimal Line-PSTH obtained from the same setof spike sequences, demonstrating the superiority of the Line-PSTH to the
Bar-PSTH in the sense of the MISE. In addition, the Line-PSTH is suitablefor comparing multiple time-dependent rates, as is the case for filtering
methods. The extrapolation method for a Line-PSTH is summarized asalgorithm 4 in the appendix.
With the extrapolation method for a Bar-PSTH (see algorithm 2), one
can also estimate how much the optimal bin size decreases (i.e., the resolution increases) with the number of spike sequences. Figure 5A shows
log-log plots of the optimal bin sizes ∗m versus m with respect to two rate
processes that fluctuate either smoothly according to the stochastic process characterized by the mean spike rate µ and the correlation of the rate

1514

H. Shimazaki and S. Shinomoto

B 60

Underlying rate

A
0

Cost function for the Bar-PSTH

0

1

2

3

0

1

2

3

2

3

C
100

Empirical cost function
Theoretical cost function

0

D

120
^

0
60

-100

^

0

0.1

0.2

0.3

0.4

0.5

0
60
^

0

0

1
Time, tFigure 2: Construction of the optimal Bar-PSTH from spike sequences. (A) Dots:
An empirical cost function, Cn (), computed from spike data according to algorithm 1. Here, 50 sequences of spikes for an interval of 20 [s] are derived
from a time-dependent Poisson process characterized by the mean rate µ and2 2
the correlation of the rate fluctuation φ(t) = σ 2 e −t /τ , with µ = 30 [spikes/s],
σ = 10 [spikes/s], and τ = 0.1 [s]. Solid line: The theoretical cost function
computed directly from the underlying fluctuating rate using equation 2.22.
(B) The underlying fluctuating rate λt . (C) Spike sequences derived from the rate.
(D) Time histograms made using three types of bin sizes: too small, optimal,
and too large.

process φ (t) = σ 2 e −t /τ or in a zigzag pattern according to the OrnsteinUhlenbeck process characterized by the mean rate µ and the correlationof the rate process φ (t) = σ 2 e −|t|/τ . These plots exhibit power law scaling
relations with distinctly different scaling exponents. The estimated exponents (−0.34 ± 0.04 for the smooth rate process and −0.56 ± 0.04 for the
zigzag rate process) are close to the exponents of m−1/3 and m−1/2 that wereobtained analytically as in equations 2.28 and 2.29 (see also Koyama &Shinomoto, 2004). In this way, we can estimate the degree of smoothness of
the underlying rate from a reasonable amount of spike data.
With the extrapolation method for a Line-PSTH (see algorithm 4 in the
appendix), the scaling relations for a Line-PSTH can be examined in a22The Bin Size of a Time Histogram

1515

A
m=30

2
1

C
Estimated critical number

^

nc

m=40

0
-1

m=60
0

4

60
8

B

40
6
20
4
2

0
0

0
1/200

1/40
1/m

1/20

20
40
60
Number of sequences, n

80

Figure 3: Extrapolation method. (A) Extrapolated cost functions Cm (|n) of
a Bar-PSTH for several values of m, computed from 30 sequences of spikes
derived from a Poisson process characterized by the mean rate µ and the cor2 2
relation of the rate fluctuation φ(t) = σ 2 e −t /τ , with µ = 30, σ = 2, and τ = 0.1.
The modestly fluctuating rate with σ = 2 is used as compared with the ratewith σ = 10 used in Figure 1. (B) The inverse optimal bin size 1/∗m plottedagainst 1/m. A solid line is a linear function fitted to the data. (C) The criticalnumber of sequences n̂c extrapolated from a smaller or larger number of spikesequences. The horizontal axis represents the number of spike sequences n used
to obtain the extrapolated cost function Cm (|n). The vertical axis represents an
extrapolated critical number n̂c . The dashed lines represent a theoretical value
nc computed using equation 2.25 and a diagonal.

similar manner. Figure 5B represents the optimal bin sizes computed for
two rate processes that fluctuate either smoothly or in a zigzag pattern.
The estimated exponents are −0.24 ± 0.04 for the smooth rate process and−0.50 ± 0.05 for the zigzag rate process. The exponents obtained by the
extrapolation method are similar to the analytically obtained exponents,
m−1/5 and m−1/2 respectively (see equations A.12 and A.13). Note that forthe smoothly fluctuating rate process, the scaling relation for the Line-PSTH
is m−1/5 , whereas the scaling relation for a Bar-PSTH is m−1/3 . In contrast,
for the rate process that fluctuates jaggedly, the exponents of the scaling
relations for both a Bar-PSTH and a Line-PSTH are m−1/2 .

1516

A

H. Shimazaki and S. Shinomoto

B

Bar-PSTH^

Line-PSTH

^

40

40

20

20

0

0

1

2

0

3

0

1

2

3

Time, tTime, tFigure 4: Comparison of the optimal Bar-PSTH and the optimal Line-PSTH.
The spike data used are the same as the data used in Figure 2. (A) The optimal
Bar-PSTH constructed from the spike sequences derived from an underlying
fluctuating rate process (dashed line). (B) The optimal Line-PSTH constructedfrom the same set of spike sequences.

A

B

Scaling relations for the Bar-PSTH
m

10

10

-1

10

-2

10

Scaling relations for the Line-PSTH
m

0

1

2

10
10
10
Number of sequences, m

3

10

-1

-2

10

0

1

2

10
10
10
Number of sequences, m

3

Figure 5: Scaling relations for the Bar-PSTH and the Line-PSTH. The optimal bin
sizes ∗m are estimated from the extrapolated cost function Cm (|n), computedfrom 100 sequences of spikes. The optimal bin sizes ∗m versus m is shown in
log-log scale. Two types of rate processes were examined: one that fluctuatessmoothly, which is characterized by the correlation of the rate fluctuation φ (t) =2 2
σ 2 e −t /τ , and the other that fluctuates jaggedly, which is characterized by thecorrelation of the rate fluctuation φ (t) = σ 2 e −|t|/τ . The parameter values are thesame as the ones used in Figure 2. (A) Log-log plots of the optimal bin sizewith respect to m for a Bar-PSTH. Lines are fitted to the data in an intervalof m ∈ [50, 500], whose regression coefficients, −0.34 ± 0.04 and −0.56 ± 0.04,
correspond to the scaling exponents for a Bar-PSTH. (B) Log-log plots of the
optimal bin size with respect to m for a Line-PSTH. Lines are fitted to the data
in an interval of m ∈ [50, 500], whose regression coefficients, −0.24 ± 0.04 and−0.50 ± 0.05, correspond to the scaling exponents for a Line-PSTH.

The Bin Size of a Time Histogram

1517

3.2 Application to Experimental Data. We also applied our method
for optimizing a Bar-PSTH to publicly available neuronal spike data (Britten, Shadlen, Newsome, & Movshon, 2004). The details of the experimental
methods are described in Newsome, Britten, and Movshon (1989) and Britten, Shadlen, Newsome, & Movshon (1992). It was reported that the neurons
in area MT exhibit highly reproducible temporal responses to identical visual stimuli (Bair & Koch, 1996). The estimation of a time-dependent rateby a PSTH is, however, sensitive to the choice of the bin size. Therefore it
would be preferable for such an appraisal of the reproducibility to be tested
with our objective method of optimizing the bin size. To make a reliable
estimate of the optimal bin size from the spike sequences of a short observation period, we took an average of the cost functions computed under
different partitioning positions.
We examine here the data recorded from a neuron under the repeatedapplication of a random dot visual stimulus with 3.2% of coherent motion(w052, nsa2004.1, Britten et al., 2004). The method for a Bar-PSTH wasapplied to the data, with close attention paid to how the optimal bin sizechanges with the number of spike sequences sampled. Figure 6 depicts theresults for the first n = 5, the first n = 20 and the total n = 50 sequences of
spikes.
The optimal bin size practically diverges (∗ = 1000 [ms]) for the n = 5
sequences, implying that with this small amount of data, a discussion about
the time-dependent response does not make sense, as far as we rely on thehistogram method. Even in this stage, it is possible to extrapolate the costfunction by means of algorithm 2. The critical number of trials, above which
the optimal bin size is finite, was estimated as n̂c ≈ 12. By increasing thenumber of spikes to n = 20, we obtained the optimal bin size ∗ = 33 [ms],implying that a discussion about the neuronal response is possible basedon this amount of data. The critical number of trials estimated from thisdata set (n = 20) is n̂c ≈ 10. The optimal bin size for the total 50 sequencesis ∗ = 23 [ms]. The critical number of sequences estimated from the total
50 sequences is n̂c ≈ 12.
Algorithm 1 confirms that the temporal rate modulation can be discussedwith a sufficiently large number of sequences, such as n = 20 or n = 50. Withalgorithm 2, three sets of sequences with n = 5, 20, and 50 suggest the criticalnumber of sequence to be n̂c ≈ 10–12.
4 Discussion
We have developed a method for selecting the bin size, so that the BarPSTH (see section 2) or the Line-PSTH (see the appendix) best represents
the (unknown) underlying spike rate. The suitability of this method wasdemonstrated by applying it to not only model spike sequences generated by time-dependent Poisson processes, but also real spike sequences
recorded from cortical area MT of a monkey.

1518

H. Shimazaki and S. Shinomoto

A
n=5

80

40
^

40

0
0

0

0.5

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4 0.6
Time, t

0.8

1

1

B
n=20
40

0
^

-40
0

0
0.05

0.1

C
n=50
0

40
^

-40
0

0.05

0.1

0

Figure 6: The Bar-PSTHs for the spike sequences recorded from a MT neuron(w052 in nsa2004.1; Britten et al., 2004). Left: The cost function (averaged overdifferent partitioning positions). Right: Spike sequences sampled and the optimal Bar-PSTH for the first 1 [s] of the recording period 2 [s]. (A) The first 5
spike sequences. (B) The first 20 spike sequences that include 5 sequences of A.
(C) The complete 50 spike sequences that include 20 sequences of B.

For a small number of spike sequences derived from a modestly fluctuating rate, the cost function does not have a minimum (Cn (T) < Cn ()for any  < T), or has a minimum at a bin size  that is comparable tothe observation period T, implying that the time-dependent rate cannot be
captured by means of the time histogram method. Our method can nevertheless extrapolate the cost function for any number of spike sequences andsuggest how many sequences should be added in order to obtain a meaningful time histogram with the resolution we require. The model data andreal data illustrated that the optimal bin size may diverge (∗ = O(T)), andeven under such conditions, our extrapolation method works reasonably
well.

The Bin Size of a Time Histogram1519

In this study, we adopted the mean integrated squared error (MISE) as
measuring the goodness of the fit of a PSTH to the underlying spike rate.
There were studies of the density estimation based on other standards suchas the Kullback-Leibler divergence (Hall, 1990) and the Hellinger distance(Kanazawa, 1993). To our knowledge, however, there are yet no practicalalgorithms based on these standards that optimize the bin size solely with
raw data. It is interesting to explore practical methods based on other standards to the estimation of the time-dependent rate and compare them with
the MISE criterion.
Use of regular (equal) bin size is a constraint particular to our method.
The regular histogram is suitable for stationary time-dependent rates. In
practical applications, however, the rate does not necessarily fluctuate ina stationary fashion but can change abruptly during a localized period oftime. In such a situation, a histogram with variable bin width can better fit
the underlying rate than a regular histogram can. In order to make a better
estimate of the underlying rate in the sense of the MISE, it is desirable
to develop an adaptive method that adjusts the bin size over time. It isalso desirable to develop the optimization method for the Bar-PSTH andthe Line-PSTH into a method for higher-order spline fittings that can becompared with filtering methods. Nevertheless, as long as the Bar-PSTHs or
the Line-PSTHs are used as conventional rate-estimation tools, our method
for selecting the bin size should be used for their construction.
Appendix: Optimization of the Line Graph Time Histogram
A line graph can be constructed by simply connecting top centers of adjacent bar graphs of the height ki /(n). Figure 7 schematically shows theconstruction of a line graph time histogram. For the same set of spike sequences, the optimal bin size (window size) of a Line-PSTH is, however,different from that of a Bar-PSTH. We develop here a method for selecting
the bin size for a Line-PSTH.
The expected heights of adjacent bar graphs for intervals of [−, 0] and[0, ] are
θ− ≡

1


θ+ ≡

1


 0
−

 

λt dt,
λt dt.

0The expected line graph L t in an interval of [− 2 , 2 ] is a line connecting top
centers of these bar graphs,
Lt =θ+ − θ−θ+ + θ−
+
t.
2


(A.1)

1520

H. Shimazaki and S. Shinomoto
+

A

0

B

^+

C
^

0
Figure 7: Line-PSTH. (A) An underlying spike rate, λt . The horizontal bars
indicate the original bar graphs θ − and θ + for adjacent intervals [−, 0] and[0, ]. The expected line graph L t is a line connecting the two top centersof adjacent expected bar graphs (− 2 , θ − ) and ( 2 , θ + ). (B) Sequences of spikes
derived from the underlying rate. (C) A time histogram for the sample sequencesof spikes. The empirical line graph L̂ t is a line connecting the two top centers ofadjacent empirical bar graphs (− 2 , θ̂ − ) and ( 2 , θ̂ + ).

The unbiased estimators of the original bar graphs are θ̂ − = k − /(n) andθ̂ + = k + /(n), where k − and k + are the numbers of spikes from n spikesequences that enter the intervals [−, 0] and [0, ], respectively. The empirical line graph L̂ t in an interval of [− 2 , 2 ] is a line connecting two top
centers of adjacent empirical bar graphs,

L̂ t =

θ̂ + + θ̂ −
θ̂ + − θ̂ −
+
t.
2


(A.2)

A.1 Selection of the Bin Size. The time average of the MISE (see equation 2.1) is rewritten by the average over the N segmented rates,

MISE =

1


 /2
−/2




E ( L̂ t − λt )2 dt.

(A.3)

The Bin Size of a Time Histogram

1521

The MISE can be decomposed into two parts:

MISE =

1


 /2




1
E( L̂ t − L t )2 dt +

−/2

 /2
−/2



(λt − L t )2 dt.

(A.4)

The first term of equation A.4 is the stochastic fluctuation of the empirical linear estimator L̂ t due to the stochastic point events, which can becomputed as
1


 /2



2
E(θ̂ + − θ + )2 .
E( L̂ t − L t )2 dt =
3
−/2


(A.5)

+
+
−
−
Here
we have  used the relations, E(θ̂ − θ )(θ̂ − θ ) = 0 and
E(θ̂ + − θ + )2 = E(θ̂ − − θ − )2 .
The second term of equation A.4 is the temporal fluctuation of λt aroundthe expected linear estimator L t . We expand the second term by insertingµ = θ −  = θ +  and obtain

1


 /2



1
(λt − L t )2 dt =

−/2
−

 /2
−/2

2





(λt − µ)2 dt

 /2
−/2

(λt − µ)(L t − µ)dt
1 /2 (L t − µ)2 dt
 −/2
 /2


1
=
(λt − µ)2 dt
 −/2



 
 
− 2 θ+ − µ θ0 − µ + 2 θ+ − µ θ∗
+

+

where
1
θ ≡

0

θ∗ ≡

 /2

2
2

−/2

λt dt, /2−/2tλt dt.

 1

2 +
(θ − µ)2 + (θ + − µ)(θ − − µ) ,
3
3

(A.6)

1522

H. Shimazaki and S. Shinomoto

To obtain the above equation, we have used relations, θ −  =θ +  = θ 0  = µ, θ ∗  = 0, (θ + − µ)2  = (θ − − µ)2 , (θ + − µ)(θ 0 − µ) =(θ − − µ)(θ 0 − µ), and (θ + − µ)θ ∗  = −(θ − − µ)θ ∗ .
The first term of equation A.6 represents a mean squared fluctuation of
the underlying rate and is independent of the choice of the bin size  (seeequation 2.8). We introduce the cost function by subtracting the variance of
the underlying rate from the original MISE:
1
Cn () ≡ MISE −
T
=

 T

(λt − µ)2 dt

0





2
2
E(θ̂ + − θ + ) − 2 (θ + − µ)(θ 0 − µ) − 2 (θ + − µ)θ ∗
3

1
2
2
(A.7)
+ E(θ + − µ) + (θ + − µ)(θ − − µ) .
3
3

The first term of equation A.7 can be estimated from the data, using thevariance-mean relation, equation 2.12. The last four terms of equation A.7
are the covariances of the expected rates θ − , θ + , θ 0 , and θ ∗ , which are notobservables. We can estimate them by using the covariance decomposition
rule for unbiased estimators:
(θ + − θ + )(θ p − θ p ) = E[(θ̂ + − E θ̂ + )(θ̂ p − E θ̂ p )]−E[(θ̂ + − θ + )(θ̂ p − θ p )],(A.8)

where p denotes −, +, 0, or ∗. The first term of equation A.8 can be computed
directly from the data, using the relation E θ̂ −  = E θ̂ +  = E θ̂ 0  = µ, andE θ̂ ∗  = 0. Unlike the Bar-PSTH, however, the mean-variance relation for
the Poisson statistics is not directly applicable to the second term. We suggest estimating these covariance terms from multiple sample sequences, as
in algorithm 3:
Algorithm 3: A Method for Bin Size Selection for a Line-PSTH
i. Divide the observation period T into N + 1 bins of width .
Count the number of spikes ki ( j) that enter ith bin from jth sequence.
(−)(+)
Define ki ( j) ≡ ki ( j) and ki ( j) ≡ ki+1 ( j) (i = 1, 2, . . . , N).
Divide the period [/2, T − /2] into N bins of width .
(0)Count the number of spikes ki ( j) that enter ith bin from jth se(∗)
quence. In each bin, compute ki ( j) ≡ 2
ti ( j)/, where ti ( j) is thetime of the th spike that enters the ith bin, measured from the centerof each bin.

( p)
( p)
( p)
ii. Sum up ki ( j) for all sequences: ki ≡ nj=1 ki ( j), where p ={−, +, 0, ∗}.

The Bin Size of a Time Histogram

1523

Average those spike counts with respect to all the bins: k̄ ( p) ≡
( p)
1 N
i=1 ki .
N(+)

Compute covariances of ki

c (+, p) ≡( p)

and ki :N
1   (+)( p)ki − k̄ (+) ki − k̄ ( p) .
N
i=1

(+)

( p)

Compute covariances of ki ( j) and ki ( j) with respect to sequences
and average over time:c̄(+, p)




N
n
(+)
( p)
ki
ki
1  1  (+)( p)≡ki ( j) −ki ( j) −
.
N
n−1
n
n
i=1

j=1

Finally, compute

σ (+, p) ≡c (+, p)c̄ (+, p)
−
.
2
(n)
n2

iii. Compute the cost function:

Cn () =2 k̄ (+)
2
1− 2σ (+, 0) − 2σ (+, ∗) + σ (+, +) + σ (+, −) .
3 (n)2
3
3

iv. Repeat i through iii while changing  to search for ∗ that minimizes
Cn ().
A.2 Extrapolation of the Cost Function. As in the case of the Bar-PSTH,
the cost function for any m sequences of spike trains can be extrapolatedusing the variance-mean relation for the Poisson statistics,

Cm (|n) =

2
3




1
1  +
E θ̂n + Cn () ,
−
m n

(A.9)

where Cn () is the original cost function (see equation A.7). With the original cost function Cn (), we can easily estimate a cost function for m sequences with algorithm 4:

1524

H. Shimazaki and S. Shinomoto

Algorithm 4: A Method for Extrapolating the Cost Function for aLine-PSTH
A. Construct the extrapolated cost function,
Cm (|n) =

2
3



1
1
−
m n



k̄ (+)+ Cn (),n2where Cn () is the cost function for the line-graph time histogram
computed for n sequences of spikes with the algorithm 3.
B. Search for ∗m that minimizes Cm (|n).
C. Repeat A and B while changing m, and plot 1/∗m versus 1/m tosearch for the critical value 1/m = 1/n̂c above which 1/∗m practicallyvanishes.
A.3 Scaling Relations of the Optimal Bin Sizes. We can obtain a theoretical cost function of a Line-PSTH directly from the mean spike rate µ andthe correlation φ(t) of the rate fluctuation. According to the mean-variancerelation based on the Cramér-Rao (in)equality (see equation 2.21), the costfunction, equation A.7, is given by
Cn () =   /2 
2µ2t2
2
φ (t1 − t2 ) dt1 dt2
1+− 2
3n  0 −/2    02
1+ 2
φ (t1 − t2 ) dt1 dt2 +φ (t1 − t2 ) dt1 dt2 .
3 0 032 0 −(A.10)

By expanding the correlation of the rate fluctuation,
 
1
1
1
φ (t) = φ (0) + φ  (0+ ) |t| + φ  (0) t 2 + φ  (0+ ) |t|3 + φ  (0) t 4 + O |t|5 ,2
624
we obtain
Cn () =181 
37 
2µ− φ (0) −φ (0+ )  +φ (0+ ) 3
3n
144
5760
49 
+φ (0) 4 + O(5 ).
2880(A.11)

Unlike the Bar-PSTH, the line graph successfully approximates the original
rate to first order in , and therefore the O(2 ) term in the cost function
vanishes for a Line-PSTH.

The Bin Size of a Time Histogram

1525

The optimal bin size ∗ is obtained from dCn () /d = 0. For a rate
process that fluctuates smoothly in time, the correlation function is a smoothfunction, resulting in φ  (0) = 0 and φ  (0) = 0 due to the symmetry φ(t) =φ(−t), and we obtain the scaling relation
∗ ∼



1280µ
181φ  (0) n1/5.

(A.12)The exponent −1/5 for a Line-PSTH is different from the exponent −1/3 for
a Bar-PSTH (see equation 2.28).
If the correlation of rate fluctuation has a cusp at t = 0 (φ  (0+ ) < 0), we
obtain the scaling relation,

 ∼ −
∗

96µ
37φ  (0+ )n1/2,(A.13)

by ignoring the higher-order terms. The exponent −1/2 for a Line-PSTH is
the same as the exponent for a Bar-PSTH (see equation 2.29).
Acknowledgments
We thank Shinsuke Koyama, Takeaki Shimokawa, and Kensuke Arai for
helpful discussions. We also acknowledge K. H. Britten, M. N. Shadlen,
W. T. Newsome, and J. A. Movshon, who have made their data available tothe public. This study is supported in part by Grants-in-Aid for Scientific
Research to S.S. from the Ministry of Education, Culture, Sports, Scienceand Technology of Japan (16300068, 18020015) and the 21st Century COE
Center for Diversity and Universality in Physics. H.S. is supported by the
Research Fellowship of the Japan Society for the Promotion of Science for
Young Scientists.
References
Abeles, M. (1982). Quantification, smoothing, and confidence-limits for single-units
histograms. Journal of Neuroscience Methods, 5(4), 317–325.
Adrian, E. (1928). The basis of sensation: The action of the sense organs. New York:
Norton.
Bair, W., & Koch, C. (1996). Temporal precision of spike trains in extrastriate cortexof the behaving macaque monkey. Neural Computation, 8(6), 1185–1202.
Blahut, R. (1987). Principles and practice of information theory. Reading, MA: AddisonWesley.
Britten, K. H., Shadlen, M. N., Newsome, W. T., & Movshon, J. A. (1992). The analysis
of visual-motion: A comparison of neuronal and psychophysical performance.
Journal of Neuroscience, 12(12), 4745–4765.

1526

H. Shimazaki and S. Shinomoto

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. Neural Signal Archive, nsa2004.1. Available online
at http://www.neuralsignal.org.
Brockwell, A. E., Rojas, A. L., & Kass, R. E. (2004). Recursive Bayesian decodingof motor cortical signals by particle filtering. Journal of Neurophysiology, 91(4),1899–1907.
Brown, E. N., Kass, R. E., & Mitra, P. P. (2004). Multiple neural spike train data
analysis: State-of-the-art and future challenges. Nature Neuroscience, 7(5), 456–
461.
Cover, T. M., & Thomas, J. A. (1991). Elements of information theory. New York: Wiley.
Daley, D., & Vere-Jones, D. (1988). An introduction to the theory of point processes. New
York: Springer-Verlag.
Dayan, P., & Abbott, L. (2001). Theoretical neuroscience: Computational and mathematicalmodeling of neural systems. Cambridge, MA: MIT Press.
DiMatteo, I., Genovese, C. R., & Kass, R. E. (2001). Bayesian curve-fitting with freeknot splines. Biometrika, 88(4), 1055–1071.
Gardiner, C. (1985). Handbook of stochastic methods: For physics, chemistry and the nauralsciences. Berlin: Springer-Verlag.
Gerstein, G. L., & Kiang, N. Y. S. (1960). An approach to the quantitative analysis ofelectrophysiological data from single neurons. Biophysical Journal, 1(1), 15–28.
Gerstein, G. L., & Perkel, D. H. (1969). Simultaneously recorded trains of action
potentials—analysis and functional interpretation. Science, 164(3881), 828–830.
Hall, P. (1990). Akaikes information criterion and Kullback-Leibler loss for histogram
density-estimation. Probability Theory and Related Fields, 85(4), 449–467.
Johnson, D. H. (1996). Point process models of single-neuron discharges. Journal ofComputational Neuroscience, 3(4), 275–299.
Kanazawa, Y. (1993). Hellinger distance and Akaike information criterion for thehistogram. Statistics and Probability Letters, 17(4), 293–298.
Kass, R. E., Ventura, V., & Brown, E. N. (2005). Statistical issues in the analysis ofneuronal data. Journal of Neurophysiology, 94(1), 8–25.
Kass, R. E., Ventura, V., & Cai, C. (2003). Statistical smoothing of neuronal data.
Network-Computation in Neural Systems, 14(1), 5–15.
Koyama, S., & Shinomoto, S. (2004). Histogram bin width selection for timedependent Poisson processes. Journal of Physics A—Mathematical and General,
37(29), 7255–7265.
Kubo, R., Toda, M., & Hashitsume, N. (1985). Nonequilibrium statistical mechanics.
Berlin: Springer-Verlag.
Newsome, W. T., Britten, K. H., & Movshon, J. A. (1989). Neuronal correlates of aperceptual decision. Nature, 341(6237), 52–54.
Révész, P. (1968). The laws of large numbers. New York: Academic Press.
Rudemo, M. (1982). Empirical choice of histograms and kernel density estimators.
Scandinavian Journal of Statistics, 9(2), 65–78.
Sanger, T. D. (2002). Decoding neural spike trains: Calculating the probability thata spike train and an external signal are related. Journal of Neurophysiology, 87(3),1659–1663.
Scott, D. W. (1979). Optimal and data-based histograms. Biometrika, 66(3), 605–610.

The Bin Size of a Time Histogram

1527

Shinomoto, S., Miyazaki, Y., Tamura, H., & Fujita, I. (2005). Regional and laminardifferences in in vivo firing patterns of primate 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.
Snyder, D. (1975). Random point processes. New York: Wiley.
Wiener, M. C., & Richmond, B. J. (2002). Model based decoding of spike trains.
Biosystems, 67(1–3), 295–300.

Received November 17, 2005; accepted September 27, 2006.