23-Mar-2007
PhD Thesis
Shimazaki, Hideaki
Kyoto University
博士(理学)
甲第12801号
理博第3111号
新制||理||1463(附属図書館)

Recipes for Selecting the Bin Sizeof a Histogram
Hideaki Shimazaki
Department of Physics, Graduate School of Science
Kyoto University, Kyoto 606-8502, Japan
December 11, 2006

Contents
1 Introduction
1.1 Optimal bin size of a histogram . . . . . . . . . . . . . .
1.2 Background and significance of the proposed methods . .

55
72 Methods
11
2.1 Selection of the bin size . . . . . . . . . . . . . . . . . . . 11
2.2 Extrapolation of the bin size . . . . . . . . . . . . . . . . 13
3 Theory
15
3.1 Construction of a bar-graph histogram . . . . . . . . . . 15
3.2 Selection of the bin size . . . . . . . . . . . . . . . . . . . 16
3.3 Extrapolation of the bin size . . . . . . . . . . . . . . . . 19
3.4 Theoretical cost function . . . . . . . . . . . . . . . . . . 21
3.4.1 Divergence of the optimal bin size . . . . . . . . . 22
3.4.2 Scaling of the optimal bin size . . . . . . . . . . . 25
4 Results
27
4.1 Selection of the bin size . . . . . . . . . . . . . . . . . . . 27
4.2 Extrapolation of the bin size . . . . . . . . . . . . . . . . 29
4.2.1 Divergence of the optimal bin size . . . . . . . . . 29
4.2.2 Scaling of the optimal bin size . . . . . . . . . . . 32
A A d-Dimensional Histogram34
A.1 Derivation of the cost function . . . . . . . . . . . . . . . 34
A.2 Selection of the bin size for a Poisson process . . . . . . . 37
1

B A Line-Graph Histogram38
B.1 Construction of a line-graph histogram . . . . . . . . . . 38
B.2 Selection of the bin size . . . . . . . . . . . . . . . . . . . 40
B.3 Extrapolation of the bin size . . . . . . . . . . . . . . . . 43
B.4 Theoretical cost function . . . . . . . . . . . . . . . . . . 43
C A Density Estimation with Sample Size Constraint48
C.1 A bar-graph histogram . . . . . . . . . . . . . . . . . . . 49
C.2 A line-graph histogram . . . . . . . . . . . . . . . . . . . 51
References

532

Abstract
A histogram is a series of bar-graphs whose heights represent the number of events within class intervals called bins. In most literature, thebin width that critically determines the goodness of the fit of the histogram to the underlying rate or density function has been selected byan individual author in an unsystematic manner. The thesis proves, in
the Poisson point process framework, that the bin width of a histogram
can be selected as the one that minimizes the formula (2k-v)/widthˆ2,
where k and v are the mean and (biased) variance of the number ofdata points in the bins.
The resolution of the histogram increases, or the optimal bin size
decreases, as the number of sampled data increases by repeating anexperimental trial. It is notable that the optimal bin size may diverge
if only a small number of experimental trials are available. In this case,any attempt to characterize the underlying function by a histogram
will lead to a spurious result. Given a paucity of data, the thesis also
provides a method that suggests how many more trials are needed until
the set of data can be analyzed with the required resolution.

3

Acknowledgements
The contents of this thesis were developed under the supervision ofProf. Shigeru Shinomoto. I should like to express my sincere thanks
to Prof. Shinomoto, who welcomed a complete stranger like me into
his group in the Nonlinear Dynamics Laboratory. His fair judgmenton the quality of research has helped me to tackle significant scientificproblems.
When I came to Kyoto three years ago, Shinsuke Koyama and Prof.
Shinomoto had already made important theoretical findings regarding
the scalings and phase transitions of the optimal histogram bin width.
Shinsuke’s pioneering work inspired me, and Prof. Shinomoto directed
my attention to developing an empirical method to choose the bin sizeof a histogram.
I enjoyed discussion on the statistics with Takeaki Shimokawa andRyota Kobayashi. Nanae Matsuno gave me many theoretical insights
needed to complete this thesis. I learned a lot about nonlinear science through discussions with Kensuke Arai, Yoji Kawamura, and Hiroya Nakao. I acknowledge Kensuke Arai for proofreading our English
manuscript submitted to journals. There were many memorable moments, including off-campus activities, that I shared with Kensuke andthe colleagues of the nonlinear dynamics group.
Before I entered the Department of Physics at Kyoto University, I
was a graduate student of the Department of Neuroscience at Johns
Hopkins University. I should like to thank my former supervisor Prof.
Ernst Niebur, who advised me to complete my Master’s degree at Hopkins. Without his help, I would not have been able to complete Ph.D
research here at Kyoto University.
Finally, I dedicate this thesis to my mother and father, who always
but implicitly encouraged me to become a scientist.
During the first two years at Kyoto University, I received a scholarship from Japan Student Services Organization. In the last year, I
was supported by the Research Fellowship of the Japan Society for the
Promotion of Science for Young Scientists.

4Chapter 1
Introduction
We introduce the problem for selecting the bin size of a histogram. The
background and significance of the solutions proposed in this thesis will
be described. The theory and methods are also available as publications
elsewhere(Shimazaki & Shinomoto, 2007a, 2007b).

1.1

Optimal bin size of a histogram

A histogram is a series of bar-graphs whose heights represent the number of events within intervals (see an example on page 6). These intervals, called bins or cells, are created by segmenting the data range ofinterest. Because of its simplicity, a histogram method is recommended
by any elementary textbook of Statistics, and is used in scientific literature. A line-graph histogram is frequently used as well. In many
applications, a histogram is used to make an inference on the underlying function that generated the data, such as a rate function or aprobability density function.
When one wishes to estimate the underlying function by a histogram, the choice of the bin size critically determines the quality of theinference. On the one hand, with a bin size that is too large, one cannot
represent details of the underlying function. On the other hand, witha bin size that is too small, the histogram fluctuates greatly and one
5

35
30
25
20
15
10
5
0
1.5

2.0

2.5

3.0

3.5

4.0

4.5

5.0

Figure 1.1: A frequency histogram for 272 observations on the durationfor eruptions of the Old Faithful geyser in Yellowstone National Park(in minutes). The bin width of the histogram was obtained by themethod proposed by this thesis (in page 12).

6

cannot discern the underlying function. For each set of data, there isan appropriate bin size which needs to be chosen based on the goodnessof the fit of the histogram to the underlying function.
The most frequently used criteria to measure the goodness of the
fit of a histogram to the underlying function is the mean integrated
squared error (MISE),ZMISE ≡ Ef {fˆn (x) − f (x)}2 dx.
(1.1)
Here f (x) is the underlying function, and fˆn (x) is a histogram constructed from the data of size n. Ef indicates expectation over an
ensemble of the histograms fˆn (x), given f (x). The thesis provides amethod to select a bin size that minimizes the MISE. The difficulty ofthe present problem comes from the fact that the underlying function
f (x) is not known.

1.2

Background and significance of the proposed methods

Theories on the optimal bin size of a histogram were mostly developedin the framework of the density estimation (Izenman, 1991). There areclassical theories that treat how the optimal bin width scales with thetotal number of data points n. It was proven that the optimal bin size
scales as n−1/3 with regard to the bar-graph density estimator (Révész,
1968; Scott, 1979). By using a Gaussian distribution as a reference,
Scott suggested to use the bin width of ∆∗ = 3.49 σ n−1/3 , where σ is astandard deviation of the samples. This approach is generally referred
to as a plug-in method. Some other extensions of Scott’s method can befound in references (Freedman & Diaconis, 1981; Wand, 1997). There
are studies on an asymptotically optimal bin size of a histogram basedon other standards. Devroye and Gyröfi studied the mean integrated
absolute error (Devroye & Györfi, 1985). Kanazawa’s approach is basedon the Hellinger distance (Kanazawa, 1993). Yet other criteria are the7Akaike’s Information Criterion (Taylor, 1987) and the Kullback-Leibler
divergence (Hall, 1990).
In the course of my graduate study, I became aware that an empirical method for choosing a histogram bin width was first presented by
Rudemo in 1982 (Eq.2.8 in Scandinavian Journal of Statistics 9: 6578). His formula is not well known by rest of the statistical community,
and is often ignored in the review papers on the histogram bin width
selection. I will demonstrate in Appendix C.1 that, if the total numberof data points sampled by experiments is predetermined, our generalformula reproduces his result. I also constructed the method of bin size
selection for a line-graph histogram under this condition, which will be
described in Appendix C.2.
All the previous studies assume that the sample size of a histogram
can be strictly controlled: a preset number of samples are randomly
picked up from a population. In many practical situations, however, thesample size is not precisely controlled by an experimentalist. Rather,indeterminate number of samples is obtained under a carefully designed
experimental protocol with a predetermined time/space constraint. A
predominant, but not exclusive, example is the data of the timing of
event occurrence in a given observation period. It is thus preferable
to have a theory that accords with these experimental designs thatfrequently appear in practice.
In this thesis, we provide a method for selecting a histogram bin
width for the data obtained through such experimental designs. The
method is applicable to selecting a size of a cell of a 2-dimensinal histogram as well. In addition to a bar-graph histogram, we also developed
a method for selecting the bin size of a line-graph histogram, which issuperior to a bar-graph in the goodness of the fit to the underlying
function.
When only a small number of experimental trials are available to
make a histogram, the estimated optimal bin size may become comparable to the range of the sampled data, implying that by constructing
a histogram, it is likely that one obtains spurious results for the estimation (Koyama & Shinomoto, 2004). Because a shortage of data
underlies this divergence, one can carry out more experiments to ob8

tain a reliable estimation. The thesis also provides the method thatcan suggest how many more experimental trials should be performed
in order to obtain a meaningful histogram with the required accuracy.
It was pointed out by Koyama and Shinomoto (Koyama & Shinomoto, 2004), and proved in general in Chapter 3, that the optimal
bin size scales differently depending on the smoothness of the underlying function. As an application of the proposed method, we also
show that the scaling relations of the optimal bin size that appears fora 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. The naive assumptions
made on the underlying function should not be accepted uncriticallybecause it is often possible to examine them from the data.
The organization of the thesis is as follows.
In Chapter 2, we provide user-friendly recipes i) for selecting thebin size of a histogram and ii) for estimating the minimum number ofexperimental trials required for constructing a histogram.
In Chapter 3, theories behind i) the method of the bin size selection
and ii) the method for estimating the minimum number of experimental
trials required for constructing a histogram will be described. We also
derive the theoretical values of i) the minimum number of experimental trials required for constructing a histogram and ii) the asymptotic
optimal bin size for a large number of experimental trials.
In Chapter 4, the empirical methods for i) selecting the bin size,ii) estimating the number of experimental trials required for histogram
construction, and iii) estimating the scaling exponents of the optimal
bin size are tested by comparing the results with those obtained by thetheoretical analysis.
In Appendix A, we generalize the theory developed in Chapter 3.
The method for selecting the bin size developed in this appendix is
applicable to a histogram with a dimension larger than one. In addition, it is applicable to selecting the bin size of a histogram whichis constructed to estimate any function related to the sampled point
events.
In Appendix B, we provide a theory and a recipe for creating a9

line-graph histogram. In addition, a theory and a recipe to estimate
the minimum number of experimental trials needed to construct a linegraph histogram will be described.
In Appendix C, the method for selecting the bin size of a bar-graphand a line-graph histogram constructed from a fixed number of samples
is described.

10Chapter 2
Methods
In this chapter, we provide a user-friendly recipe for selecting the bin
size of a bar-graph histogram of the data obtained by repeating an
identical experimental trial. We also provide the method to estimate
the minimum number of experimental trials required to construct abar-graph histogram. Recipes for a line-graph histogram are presented
in Appendix B.2.

2.1

Selection of the bin size

A bar-graph histogram is constructed simply by counting the numberof events that belong to each bin. We divide the data range of interest
into N intervals with each length given by ∆. The number of events
accumulated from all n repeated trials in the ith interval is counted aski .
The method to obtain an optimal bin size of a histogram is summarized as ‘Algorithm 1’ on page 12. In this method, the cost function
is computed by using the mean and variance of event counts within
the bins with width ∆. The bin size that minimizes this cost function
should be selected.
In the same way as in a one-dimensional histogram, a bin size ofa d-dimensional histogram can be selected by following Algorithm 1.
11

Algorithm 1: Bin size selection for a bar-graph histogram(i)

Divide the data range of interest into N bins of size ∆†, andcount the number of events ki from all n repeated trials thatenter the ith bin.

(ii)

Calculate the mean and variance of the number of events {ki }
as,
N
N
1 X
1 X
ki , and v ≡(ki − k)2 .
k≡
N i=1
N i=1

(iii) Compute the cost function,
Cn (∆) =(iv)

2k − v.
(n∆)2

Repeat i through iii while changing the bin size ∆. Find ∆∗that minimizes Cn (∆).

† For a histogram with a dimension larger than one, ∆ is a volume ofa bin.

12

In this case, the bin size ∆ is a volume of a d-dimensional cube. The
mean and variance are computed from the event counts in all the binsthat segment a total space. A d-dimensional histogram will be treatedin Appendix A.
We also constructed the method for selecting the bin size of a linegraph histogram, which is summarized as ‘Algorithm 3’ in Appendix
B.2.

2.2

Extrapolation of the bin size

The number of experimental trials plays an important role in determining the resolution of the histogram. We propose the method which
can suggest how many more experimental trials should be performed
in order to construct a meaningful bar-graph histogram with the resolution we deem sufficient. The method is summarized in Algorithm 2
on page 2.2. In this method, the cost function (in Algorithm 1) computed to obtain the bin size of a histogram for the currently available
number of experimental trials is modified to construct a cost function
for any number of experimental trials. The bin size that minimizes this‘extrapolated’ cost function is the estimated optimal bin size for thedesired number of experimental trials.
For a very small number of data, the optimal bin size may becomecomparable to the observed data range. In this case, we propose toconstruct the extrapolated cost functions for a larger number of experimental trials than the current number of trials. By examining thebin sizes that minimize these extrapolated cost functions, one can findthe number of experimental trials above which the estimated bin sizes
become smaller than the observed data range. In this way, the method
summarized in Algorithm 2 could tell how many more experimental
trials should be performed to obtain a meaningful histogram.
The extrapolation method for a line-graph histogram is summarizedas ‘Algorithm 3’ in Appendix B.3.

13

Algorithm 2: Extrapolation method for a bar-graph histogram
A Construct the extrapolated cost function,
¶
µ
1
k
1
Cm (∆|n) =−+ Cn (∆),
m n n∆2
using the mean k and variance v of the number of events obtained
from n trials. Cn (∆) is the cost function computed for n trials
with the Method 1.
B

Search for ∆∗m that minimizes Cm (∆|n). ∆∗m is the expected
optimal bin size for m trials.

C

Repeat A and B while changing m, and plot 1/∆∗m vs 1/m tosearch for the critical value 1/m = 1/n̂c above which 1/∆∗m practically vanishes.

14

Chapter 3
TheoryIn this chapter, we develop a theory on the optimal bin size of a histogram constructed from the data sampled from a Poisson point process. If the data obtained by an experimental trial is mutually independent, the sampling procedure is modeled by the inhomogeneous Poisson point process, which is the simplest point process in that the event
occurrence depends only on the instantaneous rate. Even if the dataobtained by an experimental trial is mutually dependent, as long as the
samples are obtained by repeating an identical trial, the accumulated
data are virtually independent, and obeys the inhomogeneous Poisson
point process due to a general limit theorem (Cox, 1962; Snyder, 1975;
Daley & Vere-Jones, 1988).

3.1

Construction of a bar-graph histogram

A bar-graph histogram is constructed simply by counting the numberof events that belong to each bin of width ∆. We divide the data range
T into N intervals with each length given by ∆ = T /N . The numberof events accumulated from all n sequences in the ith interval is counted
as ki . The bar height at the ith bin is given as ki /n∆. Figure 3.1 shows
the schematic diagram for the construction of a bar-graph histogram.
Given a bin of width ∆, the expected height of a bar graph for15t ∈ [0, ∆] is the time-averaged rate,Z1 ∆θ=λt dt.
∆ 0(3.1)

The total number of events k from n sequences that enter a bin of width∆ obeys a Poisson distribution with the expected number n∆θ,p(k | n∆θ) =(n∆θ)k −n∆θ
e
.
k!

(3.2)

The unbiased estimator for θ is given as θ̂ = k/(n∆), which is theempirical height of the bar graph for t ∈ [0, ∆].

3.2

Selection of the bin size

In this thesis, we assess the goodness of the fit of the estimator λ̂t to
the underlying rate λt over the data range T by the mean integrated
squared error (MISE),Z
1 T
MISE ≡E (λ̂t − λt )2 dt,
(3.3)
T 0
where E refers to the expectation over different realization of point
events, given λt .
By segmenting the range T into N intervals of size ∆, the MISE
defined in Eq.(3.3) can be rewritten asZNo2
1 ∆ 1 X n
MISE =E θ̂ (i) −λt+(i−1)∆ dt,(3.4)∆ 0 N i=1
where θ̂ (i) is an empirical height of a bar-graph histogram at ith bin.
Hereafter we denote the average over those segmented rate λt+(i−1)∆ asan average over an ensemble of (segmented) rate functions {λt } definedin an interval of t ∈ [0, ∆]:
Z
E
1 ∆D
2
(3.5)
E ( θ̂ − λt ) dt.
MISE =
∆ 0
16

A

0

B

^

C
0

Figure 3.1: The bar–graph histogram. 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 that entered each bin,
divided by the number of sequences n and the bin size ∆.
The expectation E now refers to the average over the event count,or θ̂ = k/(n∆), given a rate function λt , or its mean value, θ. The
MISE can be decomposed into two parts,
Z
E
1 ∆D
MISE =E ( θ̂ − θ + θ − λt )2 dt∆ 0D
E 1 Z ∆­®
2
= E(θ̂ − θ) +(λt − θ)2 dt.
(3.6)
∆ 0
The first and second terms are respectively the stochastic fluctuationof the estimator θ̂ around the expected mean rate θ, and the temporal fluctuation of λt around its mean θ over an interval of length ∆,17

averaged over the segments.
The second term of Eq.(3.5) can further be decomposed into twoparts,Z®
­®
1 ∆­(λt − hθi)2 dt − (θ − hθi)2 .
(3.7)
∆ 0
The first term in the rhs of Eq.(3.7) represents a mean squared fluctuation of the underlying rate λt from the mean rate hθi, and is independentof the bin size ∆, because
Z
Z
1 ∆­
1 T
2®
(λt − hθi) dt =(λt − hθi)2 dt.
(3.8)
∆ 0
T 0
We define a cost function by subtracting this term from the original
MISE,
Z
1 T
0
(λt − hθi)2 dt
Cn (∆) ≡ MISE −
T 0
D
E ­
®
2
= E(θ̂ − θ) − (θ − hθi)2 .
(3.9)
The second term in Eq.(3.9) represents the temporal fluctuation of
the expected mean rate θ for individual intervals of width ∆. As the
expected mean rate is not an observable quantity, we must replace
the fluctuation of the expected mean rate with that of the observableestimator θ̂. Using the decomposition rule for an unbiased estimator(E θ̂ = θ),DE DEE(θ̂ − hE θ̂i)2 = E(θ̂ − θ + θ − hθi)2DE ­®
= E(θ̂ − θ)2 + (θ − hθi)2 ,(3.10)
the cost function is transformed intoDE
DE
Cn (∆) = 2 E(θ̂ − θ)2 − E (θ̂ − hE θ̂i)2 .

(3.11)The second term can be taken apart into two parts,
E
E
D
D
2
2
E(θ̂ − hE θ̂i) = E (θ̂ − hθ̂i + hθ̂i − hE θ̂i)
DE
= E (θ̂ − hθ̂i)2 + E(hθ̂i − hE θ̂i)2(3.12)

18

The second term of this equation is not dependent on the choice of the0
bin size. Therefore, we add this term to the cost function Cn (∆), andredefine the cost function as0

Cn (∆) = Cn (∆) + E(hθ̂i − hE θ̂i)2
D
E
D
E
2
2
= 2 E(θ̂ − θ) − E (θ̂ − hθ̂i) .

(3.13)

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

1
E θ̂.
n∆

(3.14)

By incorporating Eq.(3.14) into Eq.(3.13), the cost function is given as
a function of the estimator θ̂,
D E
D
E
2
Cn (∆) =E θ̂ − E (θ̂ − hθ̂i)2 .
(3.15)n∆
Therefore, the cost function can be estimated by using the mean andvariance of a sample histogram,
E
2 D E D2
Cn (∆) =θ̂ − (θ̂ − hθ̂i) .
(3.16)n∆
The optimal bin size is obtained by minimizing the cost function Cn (∆):∆∗ ≡ arg min Cn (∆).
∆

(3.17)

By replacing θ̂ in Eq.(3.16) with the sample event counts, the method
is converted into a user-friendly recipe summarized in Algorithm 1.

3.3

Extrapolation of the bin size

With the method developed in the preceding subsection, we can determine the optimal bin size for a given set of experimental data. In
19this section, we develop a method to estimate how the optimal bin size
decreases when more experimental trials are added to the data set.
Assume that we are in possession of n event sequences. The fluctuation of the expected mean rate h(θ − hθi)2 i in Eq.(3.9) is replaced withthe empirical fluctuation of the histogram θ̂n using the decomposition
rule for the unbiased estimator θ̂n satisfying E θ̂n = θ,DE DEE(θ̂n − hE θ̂n i)2 = E(θ̂n − θ + θ − hθi)2
E ­D®
= E(θ̂n − θ)2 + (θ − hθi)2 .
0
The expected cost function Cm(∆) for m sequences can be obtained bysubstituting the above equation into Eq.(3.9), yielding
D
E D
E D
E
0
2
2
2
Cm (∆) = E(θ̂m − θ) + E(θ̂n − θ) − E(θ̂n − hE θ̂n i)
D
E D
E= E(θ̂m − θ)2 − E(θ̂n − θ)2 + Cn0 (∆) .
(3.18)

where Cn (∆) is the original cost function, Eq.(3.16), computed using the estimators θ̂n . We then obtain the extrapolated cost function
Cm (∆) given n sequences asD
E DECm (∆|n) = E(θ̂m − θ)2 − E(θ̂n − θ)2 + Cn (∆) .
(3.19)
Here the bias term in Eq. 3.12, which is not related to the choice of binsize, was eliminated. Using the variance-mean relation for the Poisson
distribution, Eq.(3.14), andE(θ̂m − θ)2 =we obtainµCm (∆|n) =1
1
E θ̂m =E θ̂n ,
m∆
m∆

1
1
−
m n

¶

E
1 D
E θ̂n + Cn (∆) .
∆(3.20)(3.21)

By replacing the expectation with sample event count averages, thecost function for m sequences can be extrapolated as Cm (∆|n) with20

this formula, using the sample mean k and variance v of the numbersof events, given n sequences and the bin size ∆. The extrapolation
method is summarized in Table 2.
It may come to pass that the original cost function Cn (∆) computed
for n event sequences does not have a minimum, or have a minimum
at a bin size comparable to the observed data range T . In such acase, with the method summarized in Table 2, one may estimate thecritical number of sequences nc above which the cost function has afinite bin size ∆∗ , and consider carrying out more experiments to obtaina reasonable rate estimation. In the case that the optimal bin size
exhibits continuous divergence, the cost function can be expanded as
µ
¶
1
1 1
1
Cn (∆) ∼ µ−+ u 2,(3.22)
n nc ∆∆
where we have introduced nc and u, which are independent of n. The
optimal bin size undergoes a phase transition from the vanishing 1/∆∗for n < nc to a finite 1/∆∗ for n > nc . In this case, the inverse optimal
bin size is expanded in the vicinity of nc as 1/∆∗ ∝ (1/n − 1/nc ). We
can estimate the critical value n̂c by applying this asymptotic relationˆ ∗m estimated from Cm (∆|n) for various values of m:to the set of ∆
µ
¶
1
1
1
∝
−
.
(3.23)∆∗mm n̂c
It should be noted that there are cases that the optimal bin size exhibits discontinuous divergence from a finite value. Even in such cases,
the plot of {1/m, 1/∆∗ } could be useful in exploring a discontinuoustransition from nonvanishing values of 1/∆∗ to practically vanishingvalues.

3.4

Theoretical cost function

In this section, we obtain a “theoretical” cost function directly from aprocess with a known underlying rate, λt . Note that this theoretical
21

cost function is not available in real experimental conditions in which
the underlying rate is not known.
The present estimator θ̂ ≡ k/(n∆) is a uniformly minimum varianceunbiased estimator (UMVUE) of θ, which achieves the lower bound of
the Cramér-Rao inequality (Blahut, 1987; Cover & Thomas, 1991),
" ∞
#−1
2
X
θ
∂
log
p
(k
|
n∆θ)
=
E(θ̂ − θ)2 = −p (k | n∆θ)
. (3.24)
2
∂θ
n∆
k=0
Inserting this into Eq. 3.9, the cost function is represented as®
hθi ­− (θ − hθi)2n∆Z ∆Z ∆µ1=−
φ (t1 − t2 ) dt1 dt2 ,n∆ ∆2 0 0

Cn (∆) =(3.25)

where µ = hθi is the mean rate, and φ(t1 − t2 ) = h(λt1 − µ)(λt2 − µ)i is
the autocorrelation function of the rate fluctuation, λt − µ.

3.4.1

Divergence of the optimal bin size

We obtain analytical solution of the critical number nc , below which
the optimal bin size that minimizes the cost function diverges. The
cost function with a large bin size can be rewritten asZ ∆1µ(∆ − |t|)φ(t) dt−
Cn (∆) =n∆ ∆2 −∆
Z
Z ∞
µ
1 ∞
1
∼
−
φ(t) dt + 2|t|φ(t) dt,(3.26)n∆ ∆ −∞∆ −∞based on the symmetry φ(t) = φ(−t) for a stationary process. Eq. 3.26
can be identified with Eq. 3.22 with parametersÁZ ∞φ(t) dt ,nc = µ(3.27)−∞Z ∞u=|t|φ(t) dt.
(3.28)−∞

22

The minimum number of repeated trials required for the constructionof a histogram is given by dnc e. For the process giving the correlation of
the form φ (t) = σ 2 e−|t|/τ and the process giving a Gaussian correlation
22
φ(t) = σ 2 e−t /τ √, the critical numbers are obtained as nc = µ/2σ 2 τand nc = µ/σ 2 τ π, respectively.
Note that the expansion of the cost function with respect to 1/∆
assumes that the optimal bin size takes an infinitely large value at theonset of the transition. However, the finite optimal bin size may appear
suddenly as we increase the number of trials. The transitions of the
optimal bin width from infinite to finite length can be systematicallyexamined by considering Fourier domain of the cost function.
The Fourier transform of the window function, whose height is 1/∆
within the window of width ∆ and zero outside the window, is a sincfunction given byZ ∆/2
1 −jωtsin(ω∆/2)H∆ (ω) =edt =.
(3.29)
ω∆/2−∆/2 ∆
Let the Fourier transform of the correlation φ(t) be Φ(ω) (the power
spectrum of the rate process). Using H∆ (ω) and Φ(ω), the cost function
Eq. 3.25 is written as
Z ∞
n
o
1
2 µ
H∆ (ω)− Φ(ω) dω.
(3.30)
Cn (∆) =2π −∞nThe cost function vanishes for the large bin size: lim∆→∞ Cn (∆) = 0.
Therefore, if the finite optimal bin size ∆∗ exists, it should satisfy
Cn (∆∗ ) < 0.
We first examine the cost function of the stochastic processes whose
power spectrum has a maximum at origin (Φ(0) > Φ(ω) for ω 6= 0).
For a small number of sequences that satisfies µ/n > Φ(0), the costfunction is positive for all ∆. Therefore the optimal bin size diverges
if n < µ /Φ(0). The critical number of trial is then given by nc =µ /Φ(0), which is identical to Eq.3.27. To examine the condition near
the transition point, suppose that n is slightly larger than nc . The termµ/n − Φ(ω) in the integral of Eq.3.30 is negative in the vicinity of the
23

D

A

2
0.5
1
0

‑5

0

0

5

‑5

0

5

E

B
10

10

5

5

0

0

‑5

0

5

10

15

‑5

20

0

5

10

15

20

F

C

8
2

6
4

1

2
0
0.02

0.04

0.06

0.08

0.1

0
0.02

0.04

0.06

0.08

1/n

1/n

Figure 3.2: The first and the second order phase transitions observed for
optimal bin sizes with two different rate processes. (Left; A, B, C) a pro2
2
cess whose mean is µ = 30 and correlation function is φ(t) = σ 2 e−t /τ(σ 2 = 9, τ = 0.1) (Right; D, E, F) a process whose mean is µ = 30 andcorrelation function is φ(t) = σ 2 (1 − |t|/τ )e−|t|/τ (σ 2 = 25, τ = 0.1).
(Top; A, D) Power spectra of the rate process. (Middle; B, E) Theoretical cost functions derived from Eq.3.25 (B n = 10, 40,and 100; E
n = 10, 20,and 30). The asterisks indicate optimal bin sizes (Bottom,CF) The inverse of optimal bin sizes as a function of the inverse of n.
The dots represent the optimal bin sizes for integer n.

24origin and positive elsewhere. As the bin size ∆ increases, the spectraof H∆ (ω)2 vanish except for the origin whose height is kept constant as
1. Therefore the narrow negative spectrum range of µ/n − Φ(ω) can be
exploited by H∆ (ω)2 with a large bin size: An appropriately large bin
size yields the negative cost function. This indicates that the optimal
bin size appears as an infinitely large ∆. This type of phase transitions
is called the second order (continuous) phase transitions. For exam22ple, a process with a Gaussian correlation function (φ(t) = σ 2 e−t /τ )displayed in Fig. 3.4.1A exhibit the second order phase transition.
For the stochastic processes whose power spectrum has a maximum
other than origin, it is possible that the first order (discontinuous) phasetransitions take place. Figure 3.4.1A displays the power spectrum of
the process with a correlation function φ(t) = σ 2 (1 − |t|/τ )e−|t|/τ . This
process has a negative correlation for |t| > τ , and the power spectrum at
origin is not a maximum. The cost function analytically obtained from
Eq. 3.30 (Figure 3.4.1E) reveals that, when the number of sequences isincreased, two local minima appear (one at origin and the other at a
finite bin size). The optimal bin size undergoes a discontinuos phasetransition as shown in Figure 3.4.1F. The example is an extreme case,where the first order phase transition always occur because Φ(0) = 0.
It should be noted that even if the power spectrum has a maximum
other than origin, the second order phase transition could take place,
too.

3.4.2

Scaling of the optimal bin size

Based on the same theoretical cost function Eq. 3.25, we can also derivescaling relations of the optimal bin size, which is achievable with a large
number of event sequences. With a small ∆, the correlation of rate
fluctuation φ(t) can be expanded as¡ ¢1
φ (t) = φ (0) + φ0 (0+ ) |t| + φ00 (0) t2 + O |t|3 ,2and we obtain an expansion of Eq. 3.25 with respect to ∆,11µ− φ (0) − φ0 (0+ ) ∆ − φ00 (0) ∆2 + O(∆3 ). (3.31)
Cn (∆) =
n∆
3
12
25

The optimal bin size ∆∗ is obtained from dCn (∆) /d∆ = 0. Fora rate that fluctuates smoothly in time, the correlation function is asmooth function of t, resulting in φ0 (0) = 0 due to the symmetry φ(t) =φ(−t). In this case, we obtain the scaling relation
µ
∗

∆ ∼

6µ
− 00
φ (0) n¶1/3.

(3.32)For a rate that fluctuates in a zigzag pattern, in which the correlation of rate fluctuation has a cusp at t = 0 (φ0 (0+ ) < 0), we obtain thescaling relation
µ
¶1/2
3µ
∗
∆ ∼ − 0,(3.33)
φ (0+ )nby ignoring the second order term in Eq.3.31. Examples that exhibit
this type of scaling are the Ornstein-Uhlenbeck process and the random
telegraph process (Kubo & Hashitsume, 1985; Gardiner, 1985).

26

Chapter 4
Results
In this chapter, we apply the presently developed methods for the bargraph histogram and the line-graph histogram to sequences of events
generated by the simulations of inhomogeneous Poisson point processes.

4.1

Selection of the bin size

We applied the method for estimating the optimal bin size of a bargraph histogram to a set of event sequences derived from an inhomogeneous Poisson process. The rate of the Poisson point process usedin the simulation is a stationary stochastic process. Figure 4.1A showsthe ‘empirical’ cost function Cn (∆) computed from a set of n event
sequences according to Algorithm 2.1. The empirical cost function
approximates the ‘theoretical’ cost function for the rate process, according to Eq. 3.25 in Section 3. In Figure 4.1D, a histogram withthe optimal bin size is compared with those with non-optimal bin sizes,
demonstrating the effectiveness of optimizing the bin size. We also constructed a method for selecting bin size of a line-graph histogram, which
is summarized as ‘Algorithm 3’ in Appendix B.2. A line graph can beconstructed by simply connecting top-centers of adjacent bar graphs.
Figure 4.2 compares the optimal bar-graph histogram and the optimalline-graph histogram, obtained from the same set of event sequences,
27

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 4.1: Construction of the optimal bar-graph histogram from point
events. A: (Dots): an ‘empirical’ cost function, Cn (∆), computed fromthe data according to the Algorithm 1. Here, 50 sequences of events
for an interval of 20 are derived from an inhomogeneous Poisson process characterized by the mean rate µ and the correlation of the rate22
fluctuation φ(t) = σ 2 e−t /τ , with µ = 30, σ = 10, and τ = 0.1. (Solidline): the ‘theoretical’ cost function computed directly from the underlying fluctuating rate using Eq. 3.25. B: the underlying fluctuating rate
λt . C: event sequences derived from the rate. D: bar-graph histograms
made using three types of bin sizes: too small, optimal, and too large.

28

demonstrating the superiority of the line-graph histogram to the bargraph histogram in the sense of the MISE. In addition, the line-graphhistogram is suitable for the comparison of multiple rates, as is the casefor filtering methods.

4.2

Extrapolation of the bin size

The optimal bin size for the number of experimental trials m can be
estimated from the currently available data of size n by minimizing
the extrapolated cost function in the Algorithm 2. In this method,
the cost function of m sequences is obtained by modifying the original
cost function, Cn (∆), computed from the event count statistics of nsequences of events.

4.2.1

Divergence of the optimal bin size

For a small number of sampled data, the estimated optimal bin size
may diverge, implying that by constructing a histogram, it is likelythat one obtains spurious results for the estimation (Koyama & Shinomoto, 2004). Because a shortage of data underlies this divergence,
one can carry out more experiments to obtain a reliable estimation.
With the Algorithm 2, one can estimate how many more experimental
trials should be performed in order to construct a meaningful bar-graphhistogram.
We applied this extrapolation method for a bar-graph histogram toa set of event sequences derived from the smoothly regulated Poissonprocess. Figure 4.3A depicts the extrapolated cost function Cm (∆|n)
for several values of m, computed from a given set of n(= 30) event
sequences. Figure 4.3B represents the dependence 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 linearly with m for 1/m ≤1/ n̂c . By fitting the linear function, we estimated the critical numberof sequences n̂c . In Figure 4.3C, the critical number of sequences n̂cestimated from a smaller or larger n is compared to the theoretical value
29

A

A bar-graph histogram

^

40
20
0

B

0

1

2

3

A line-graph histogram

^

40
20
0

0

1

2

3

Figure 4.2: The comparison of the optimal bar-graph histogram andthe optimal line-graph histogram. The data of the point events are thesame as the data used in Fig. 4.1. A: the optimal bar-graph histogram
constructed from the event sequences derived from an underlying fluctuating rate process (dashed line). B: the optimal line-graph histogram
constructed from the same set of event sequences.

30

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 sequence s, n

80

Figure 4.3: The extrapolation method. A: extrapolated cost functions
Cm (∆|n) of a bar-graph histogram for several values of m, computedfrom 30 sequences of events derived from a Poisson process characterized by the mean rate µ and the correlation of the rate fluctuation
22
φ(t) = σ 2 e−t /τ , with µ = 30, σ = 2, and τ = 0.1. The modestly fluctuating rate with σ = 2 is used as compared with the rate with σ = 10
used in Fig. 1. B: the inverse optimal bin size 1/∆∗m plotted against
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 numberof event sequences. The horizontal axis represents the number of event
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 Eq. 3.27,
and a diagonal.

31

of nc computed from Eq. 3.27. It is noteworthy that the estimated n̂capproximates the theoretical nc well even from a fairly small numberof event sequences, with which the estimated optimal bin size diverges.

4.2.2

Scaling of the optimal bin size

With the extrapolation method for a bar-graph histogram (Algorithm2), one can also estimate how much the optimal bin size decreases (i.e.,
the resolution increases) with the number of event sequences.
Figure 4.2.2A shows log-log plots of the optimal bin sizes ∆∗m versus
m with respect to two rate processes that fluctuate either smoothly orin a zigzag pattern. 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 thatwere obtained analytically as in Eqs. 3.32 and 3.33. In this way, we
can estimate the degree of smoothness of the underlying rate from areasonable amount of data.
With the extrapolation method for a line-graph histogram (Algorithm 4 in Appendix B.3), the scaling relations for a line-graph histogram can be examined in a similar manner. Figure 4.2.2B represents
the optimal bin sizes computed for two rate processes that either fluctuate 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 aresimilar to the analytically obtained exponents, m−1/5 and m−1/2 respectively (see Eqs. B.15 and B.16). Note that for the smoothly fluctuating
rate process, the scaling relation for the line-graph histogram is m−1/5 ,whereas the scaling relation for a bar-graph histogram is m−1/3 . In
contrast, for the rate process that fluctuates jaggedly, the exponents ofthe scaling relations for both a bar-graph histogram and a line-graphhistogram are m−1/2 .

32

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 sequence s, m

3

10

-1

-2

10

0

1

2

10
10
10
Number of sequences, m

3

Figure 4.4: Scaling relations for the bar-graph histogram and the linegraph histogram. The optimal bin sizes ∆∗m are estimated from theextrapolated cost function Cm (∆|n), computed from 100 sequencesof 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 fluctua2
2tion φ (t) = σ 2 e−t /τ , and the other that fluctuates jaggedly, which ischaracterized by the correlation of the rate fluctuation φ (t) = σ 2 e−|t|/τ .
The parameter values are the same as the ones used in Fig. 4.1. A:
log-log plots of the optimal bin size with respect to m for a bar-graphhistogram. Lines are fitted to the data in an interval of m ∈ [50, 500],
whose regression coefficients, −0.34±0.04 and −0.56±0.04, correspondto the scaling exponents for a bar-graph histogram. B: log-log plots of
the optimal bin size with respect to m for a line-graph histogram. 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-graph histogram.

33

Appendix A
A d-Dimensional Histogram
In this Appendix, we generalize the method for selecting the bin sizeof a histogram developed in Chapter 3. In Appendix A.1, we re-derive
the cost function so that the formula is applicable to selecting the bin
size of a histogram of d-dimension, where d is positive integer. In addition, the cost function derived in this appendix is generally applicable
to a histogram that depicts the function related to point events. In
Appendix A.2, we explicitly derive the cost function for a histogramfor density estimation.

A.1

Derivation of the cost function

For a histogram of a multivariate datasets x, MISE is defined by multiple integral:Zno2
1
MISE =
E fˆ (x) − f (x) dx.
(A.1)|V | V
where V is the d-dimensional data space of interest and |V | is its volume. The average value of a function f (x) in the ith bin is given
by
Z
1
¯
fi =
f (x) dx.
∆i Vi
34

where Vi represents the space of ith bin, and ∆i is the volume of theith bin (∆i = |Vi |). As is in one dimensional histogram, we divided themultiple integral into the N segmented bins:
Z
o2
n
1
ˆ
MISE =
Ef f (x) − f (x) dx|V | V
Z
N
n
o2
1 X 1
=
Efi fˆi −f (x) dx.
N i=1 ∆i Vi(A.2)

By denoting the average over bins with a bracket, we rewrite MISE as
¿ Z
n
o2 À
1
MISE =
Ef fˆi −f (x) dx .
(A.3)
∆i Vi i
By inserting the ideal model f¯i , the MISE is decomposed into two parts,
¿ Z
n
o2 À
1
ˆ
¯
¯
Ef fi −f i +f i −f (x) dx
MISE =(A.4)∆i Vi iÀ¿³
´2 À ¿ 1 Z ©
ª2¯
ˆ¯f (x) −f i dx .
= Efi fi − fi
+
∆i Vi
R
By inserting the total average µ = V1 V f (x)dx, the second term canbe further decomposed into two,
¿ Z
À
À ¿ Z
À ¿ Z
ª2
©
ª2
©
1
1
1
2
¯
¯
fi −µ dxf (x) −f i dx ={f (x) − µ} dx −
∆i Vi
∆i Vi
∆i Vi
Z
D¡
¢2 E
1
.
{f (x) − µ}2 dx− f¯i −µ=|V | V(A.5)
where the first term is irrespective of the choice of the bin size. We
thus introduce the cost function,
Z
1
0
C ({∆i }) = MISE −{f (x) − µ}2 dx|V | VD
E D¡¢2 E= Efi (fˆi − f¯i )2 − f¯i −µ
.
(A.6)
35

From the decomposition rule, we have
D
D
E E D
E D
¡
­ ®¢2 EEfi (fˆi − Efi fˆi )2 = Efi (fˆi − f¯i )2 + Efi f¯i − f¯i. (A.7)
By substituting this equation, the cost function becomes
¿³
¿
³
´2 À
D
E´2 À
C ({∆i }) = 2 Efi fˆi − f¯−Effˆi − E fˆi
.

(A.8)

The second term is written as¿³¿³¿³D E DDE´2 À
D E´2 ÀE´2 ÀEffˆi − Efi fˆi= Effˆi − fˆi+Ef
fˆi − Efi fˆi.
(A.9)
The second term of the above equation does not depend on the choiceof the bin size ∆. We add this term to the cost function, and redefine
the cost function as¿³D E DE´2 À0ˆˆC ({∆i }) = C ({∆i }) + Ef
fi − Efi fi
¿
¿³
³
´2 À
D E´2 Àˆ¯
ˆ= 2 Efi fi − f−Ef
fi − fˆi.
(A.10)
Note that throughout the derivation of the cost function, we only
assumed that each height of the histogram is unbiased estimator of theaverage of the underlying function in the bin. The cost function given
by A.10 is applicable to not only obtaining the bin size of histogramfor a density (or rate) estimation, but also obtaining the bin size ofa histogram that depicts the function related to a point process, suchas an auto-correlation function of the underlying rate. When the experimental trials are repeated to obtain the data, the first term can beevaluated as the variance of the histograms constructed from individual
trials.

36

A.2

Selection of the bin size for a Poissonprocess

For a rate estimation of a Poisson point process, the bar-height is defined askifˆi = θ̂i =,(A.11)
n∆i
where n is the number of repeated trials. The number of events ki inthe volume Vi follows a Poisson distribution (Diggle, 1983) with mean,Znf (x) dx =n∆fi .
(A.12)
Vi

Since the mean and variance of the Poisson point process equals, we
have
¿ ³
À
´2 À ¿ 1
E fˆi .
E fˆi − f¯i(A.13)=n∆i
The cost function is then written as¿À¿³
D E´2 À1 ˆ
2
ˆCn ({∆i }) = Ef
fi −Effi − fˆi.
n∆i(A.14)

By replacing the expectation by one sample histogram, we obtain*µ¿¿ À¶2 +
À
1
ki
ki
1
2
k − 2
−
.
Cn ({∆i }) = 2
2 i
n
∆i
n
∆i
∆i
For a regular histogram, whose volume is equal in all bins, we have
Cn (∆) =

1
n2 ∆2

©

­
®ª
2 hki i − (ki − hki i)2 .

37

Appendix B
A Line-Graph Histogram
In this Appendix, we provide a thoery and a recipe for i) selecting thebin size of a line-graph histogram and ii) estimating required numberof experimental trials to construct a line-graph histogram.

B.1

Construction of a line-graph histogram

A line graph can be constructed by simply connecting top-centers ofadjacent bar graphs of the height ki /(n∆). Figure B.1 schematicallyshows the construction of a line-graph histogram. For the same setof event sequences, the optimal bin size (window size) of a line-graphhistogram is, however, different from that of a bar-graph histogram.
We develop here a method for selecting the bin size for a line-graphhistogram.
The expected heights of adjacent bar graphs for intervals of [−∆, 0]
and [0, ∆] are
Z
1 0
−
θ ≡
λt dt,
∆ −∆
Z
1 ∆
+
θ ≡
λt dt.
∆ 0

38

+

A

0

B

^+

C
^

0

Figure B.1: The line-graph histogram. A: an underlying rate, λt . The
horizontal bars indicate the original bar graphs θ− and θ+ for adjacentintervals [−∆, 0] and [0, ∆]. The expected line graph Lt is a line connecting the two top-centers of adjacent expected bar graphs (− ∆2 , θ− )and ( ∆2 , θ+ ). B: sequences of events derived from the underlying rate.
C: a line-graph histogram for the sample sequences of events. The empirical line graph L̂t is a line connecting the two top-centers of adjacentempirical bar graphs (− ∆2 , θ̂− ) and ( ∆2 , θ̂+ ).

39

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

(B.1)

The unbiased estimators of the original bar graphs are θ̂− = k − /(n∆)
and θ̂+ = k + /(n∆), where k − and k + are the numbers of events from n
event sequences 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 =

B.2

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

(B.2)

Selection of the bin size

The time average of the MISE (Eq. 3.3) is rewritten by the averageover the N segmented rates,
ZE
1 ∆/2 D2
E (L̂t − λt ) dt.
(B.3)
MISE =∆ −∆/2The MISE can be decomposed into two parts,
ZZ
E®
1 ∆/2 D1 ∆/2 ­
2
MISE =
E(L̂t − Lt ) dt +(λt − Lt )2 dt. (B.4)∆ −∆/2∆ −∆/2
The first term of Eq. B.4 is the stochastic fluctuation of the empiricallinear estimator L̂t due to the stochastic point events, which can becomputed asZE
1 ∆/2 D
2
E(L̂t − Lt ) dt
∆ −∆/2
E 1D
E
1D
(B.5)
E(θ̂+ − θ+ )2 +E(θ̂− − θ− )2 .
=
3
3
40

Here we have used the relations, E(θ̂+ − θ+ )(θ̂− − θ− ) = 0.
The second term of Eq. B.4 is the temporal fluctuation of λt aroundthe expected linear estimator Lt . We expands the second term byinserting µ = hθ− i = hθ+ i = hθ0 i, and obtain
1
∆

Z ∆/2

­

(λt − Lt )

2®

−∆/2

1
dt =
∆
−

2
∆

+

Z ∆/2

­

®
(λt −µ)2 dt(B.6)

h(λt −µ) (Lt −µ)i dt

(B.7)

−∆/2
Z ∆/2

1
∆

−∆/2
Z ∆/2
­

®
(Lt −µ)2 dt

(B.8)

−∆/2

The second term B.7 is computed as
­¡
¢¡
¢® ­¡
¢¡
¢®
(B.7) = θ+ −µ θ0 −µ + θ− − µ θ0 − µ­¡¢ ® ­¡¢ ®
+ θ+ −µ θ∗ − θ− −µ θ∗ ,
where
1
θ ≡
∆
0

θ∗ ≡

Z ∆/2

2
∆2

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

To obtain the above equation, we have used relations, hθ− i = hθ+ i =hθ0 i = µ. The third term B.8 can be computed as(B.8) =

¢2 E 1 ­¡ +¢¡¢® 1 D¡ − ¢2 E
1 D¡ +
θ −µ+
θ − µ θ− − µ +
θ −µ
.
3
3
3

The first term B.6 of the systematic error represents a mean squared
fluctuation of the underlying rate, and is independent of the choiceof the bin size ∆ (See Eq. 3.8). We introduce the cost function by41

subtracting the variance of the underlying rate from the original MISE:
Z
1
0
Cn (∆) ≡ MISE −(λt − µ)2 dt
T
E 1D
E
1D
E(θ̂+ − θ+ )2 +E(θ̂− − θ− )2=3­¡¢ ¡ 0 ¢® 3­¡ −¢¡¢®
+− θ −µ θ −µ − θ − µ θ0 − µ­¡¢ ® ­¡¢ ®
− θ+ −µ θ∗ + θ− −µ θ∗¢2 E 1 D ¡ − ¢2 E
1 D¡ +
+
θ −µ
+
θ −µ
3
3
¢¡
¢®
1 ­¡ ++
θ − µ θ− − µ .
(B.9)
3
D
E D
E
+
+ 2
−
− 2
By using the relations, E(θ̂ − θ ) = E(θ̂ − θ ) , h(θ+ − µ)2 i =h(θ− − µ)2 i, h(θ+ − µ) (θ0 − µ)i = h(θ− − µ) (θ0 − µ)i, and h(θ+ − µ) θ∗ i =− h(θ− − µ) θ∗ i, the cost function becomes
E
2D
+
+ 2
E(θ̂ − θ )Cn (∆) =
3­¡
¢¡
¢®
­¡
¢ ®
− 2 θ+ −µ θ0 −µ − 2 θ+ −µ θ∗¢2 E 1 ­¡ +¢¡¢®
2 D¡ ++
θ −µ+
θ − µ θ− − µ .
(B.10)
3
3
The first term of Eq. B.10 can be estimated from the data, using
the variance-mean relation (Eq. 3.14). The last four terms of Eq. B.10
are the covariances of the expected rates θ− , θ+ , θ0 , and θ∗ , which are
not observables. We can estimate them by using the covariance decomposition rule for unbiased estimators:
D ³
D
E´ ³
D
E´E
E θ̂+ − E θ̂+θ̂p − E θ̂p
D h³
´³
´iE ­
®
= E θ̂+ − θ+ θ̂p − θp
+ (θ+ − hθ+ i)(θp − hθp i) . (B.11)
where p denotes −, +, 0, or ∗. Unlike the bar-graph histogram, however, the mean-variance relation for the Poisson statistics is not directly
applicable to the second term. We suggest estimating these covariance42

terms from multiple sample sequences. The r.h.s of Eq. B.11 is decomposed into
D³
D
E´ ³
D
E´E
+
+
p
p
E θ̂ − E θ̂θ̂ − E θ̂
D³
D E´ ³D E´E³D E ­ ®´ ³D E´
= E θ̂+ − θ̂+θ̂p − θ̂p
+ E θ̂+ − θ+θ̂p − hθp i .
The second term of the r.h.s. of this equation is not dependent on thebin size selection. Therefore, we redefine the cost function appropriately
as³D E ­ ®´ ³D E ­ ®´
Cn (∆) = Cn0 (∆) + 2E θ̂+ − θ+θ̂0 − θ0
´ 2 ³D E ­ ®´2³D E ­ ®´ ³D E++∗∗
+ 2E θ̂ − θθ̂ − hθ i − E θ̂+ − θ+3This cost function can be estimated by following a recipe summarizedas ‘Algorithm 3’.

B.3

Extrapolation of the bin size

As in the case of the bar-graph histogram the cost function for any msequences of event sequences can be extrapolated using the variancemean relation for the Poisson statistics,µ¶
2
1
1 D +ECm (∆|n) =(B.12)E θ̂n + Cn (∆) ,−3∆ m nwhere Cn (∆) is the original cost function (Eq. B.10). With the original cost function Cn (∆), we can easily estimate a cost function for m
sequences with ‘Algorithm 4’.

B.4

Theoretical cost function

We can obtain a ‘theoretical’ cost function of a line-graph histogram
directly from the mean rate µ and the correlation φ(t) of the rate fluc43

tuation. According to the mean-variance relation based on the CramérRao (in)equality (Eq. 3.24), the cost function (Eq. B.10) is given by¶Z ∆ Z ∆/2 µ2µ22t2
Cn (∆) =−1+
φ (t1 − t2 ) dt1 dt2
3n∆ ∆2 0 −∆/2∆
Z ∆Z ∆Z ∆Z 0
2
1+
φ (t1 − t2 ) dt1 dt2 +φ (t1 − t2 ) dt1 dt2 .
3∆2 0 03∆2 0 −∆
(B.13)

Scaling relations of the optimal bin sizes
By expanding the correlation of the rate fluctuation,
¡ ¢
1
1
1
φ (t) = φ (0)+φ0 (0+ ) |t|+ φ00 (0) t2 + φ000 (0+ ) |t|3 + φ0000 (0) t4 +O |t|5 ,2
624
we obtain
37 0
181 000
49 0000
2µ−φ (0)−
φ (0+ ) ∆+φ (0+ ) ∆3 +φ (0) ∆4 +O(∆5 ).
Cn (∆) =
3n∆
144
5760
2880
(B.14)
Unlike the bar-graph histogram, 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-graph histogram.
The optimal bin size ∆∗ is obtained from dCn (∆) /d∆ = 0. For arate process that fluctuates smoothly in time, the correlation function
is a smooth function, resulting in φ0 (0) = 0 and φ000 (0) = 0 due to thesymmetry φ(t) = φ(−t), and we obtain the scaling relation
µ
¶1/5
1280µ
∗
∆ ∼
.
(B.15)
181φ0000 (0) n
The exponent −1/5 for a line-graph histogram is different than the
exponent −1/3 for a bar-graph histogram (Eq. 3.32).
If the correlation of rate fluctuation has a cusp at t = 0 (φ0 (0+ ) < 0),
we obtain the scaling relation
µ
¶1/2
96µ
∗
,
(B.16)∆ ∼ −37φ0 (0+ )n44

by ignoring the higher order terms. The exponent −1/2 for a linegraph histogram is the same as the exponent for a bar-graph histogram(Eq. 3.33).

45

Algorithm 3: Bin size selection for a line-graph histogram(i)

Divide the observation period T into N + 1 bins of width ∆.
Count the number of events 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 events ki (j) that enter ith bin from jth
sequence.
P(∗)
In each bin, compute ki (j) ≡ 2 ` t`i (j)/∆, where t`i (j) is thetime of the `th event that enters the ith bin, measured fromthe center of each bin.

(ii)

(p)
Sum up ki (j) for all sequences:

(p)
ki ≡

n
X

(p)

ki (j), where p =j=1{−, +, 0, ∗}.
Average those event counts with respect to all the bins: k̄ (p) ≡N
1 X (p)k .
N i=1 i(+)

Compute covariances of kic(+, p)(p)and ki :N´³´
1 X ³ (+)
(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:
Ã
!Ã
!
N
n
(+)
(p)
X
X
1
1
k
k
(+)
(p)
c̄(+, p) ≡ki (j) − iki (j) − i.
N i=1 n − 1 j=1nn

46Algorithm 3: continued from page 46.
Finally, compute
σ (+, p) ≡c̄(+, p)c(+, p)
−
.
(n∆)2
n∆2

(iii) Compute the cost function,
Cn (∆) =(iv)

2 k̄ (+)
2
1− 2σ (+, 0) − 2σ (+, ∗) + σ (+, +) + σ (+, −) .
23 (n∆)
3
3

Repeat i through iii while changing ∆. Find ∆∗ that minimizes
Cn (∆).

Algorithm 4: Extrapolation method for a line-graph histogram
A Construct the extrapolated cost function,µ¶
2 1
1 k̄ (+)
Cm (∆|n) =−+ Cn (∆),
3 m n n∆2where Cn (∆) is the cost function for the line-graph histogram
computed for n trials with the Algorithm 3.
B

Search for ∆∗m that minimizes Cm (∆|n).

C

A and B while changing m, and plot 1/∆∗m vs 1/m to search
for the critical value 1/m = 1/n̂c above which 1/∆∗m practicallyvanishes.

47

Appendix C
A Density Estimation with
Sample Size Constraint
In this Appendix, we derive cost functions of a bar- and a line-graphhistogram when the data with predetermined size is sampled from a
population. In Chapter 3, we investigated the histogram of the datathat are derived from an inhomogeneous Poisson point process, where
the sample size is not determined. Due to the assumption of iid randomvariables, histograms of the data obtained by the two different experimental designs can be equally seen as the estimation of the probability
density function. However, the difference between the two experimental designs appears when one considers ensemble of histograms as in
the MISE criterion.
The MISE is composed of the sampling and systematic errors asshown in Eq.3.6. A similar equation for a line-graph histogram wasgiven by Eq.B.4. For the Poisson point process, the sampling error simply obeys the variance of the Poisson statistics. When a fixed numberof samples are collected from a population, the sampling error becomes
smaller than the variance of the Poisson statistics. As a result, theoptimal bin size obtained under the sample size constraint is slightly
smaller than that obtained under the Poisson point process framework.

48

C.1

A bar-graph histogram

The bar height of a histogram for a density estimation is given byθ̂n = ki /n∆. In this context, n should be read as the number of samples(it is not the number of experimental trials). We label the sequence ofthe samples by j (= 1, · · · , n), and introduce ui (j). ui (j) is a unity if
the jth sampled data is observed in the ith bin (i = 1, · · · , N ) , andzero otherwise. Summation of ui (j) with respect to j equals ki :
ki =

n
X

ui (j) .

(C.1)

j=1

The sampling error for the ith bin is given by()2nn
1 1 X ui (j) 1 X ui (j)−n n − 1 j=1∆
n j=1 ∆
By using the notation, ūi = n1
n
X

n
P

(C.2)

ui (j) = ki /n, we obtain

j=1

1
1
1
(ui (j) − ūi )2 =2n∆ n − 1 j=1
n∆2

(

n
X

)

1
n
ui (j)2 −ū2in − 1 j=1
n−1
µ
¶
1
1
1
2
=
ki −
k
(C.3)
n∆2 n − 1n (n − 1) i

To obtain the last equality, we have used a relation,
ki =

n
X

ui (j) =j=1nXui (j)2 .

(C.4)

j=1

The last equality holds because ui (j) = {0, 1}. This relation manifests
the constraint on the total sample size.
In the derivation of the cost function Eq.3.13 in Section 3.2 in Chapter 3, we have only assumed that a histogram is an unbiased estimator49

of the underlying function. Therefore, the cost function in the form of
Eq.3.13,DE
DE
Cn (∆) = 2 E(θ̂ − θ)2 − E (θ̂ − hθ̂i)2 .
(3.13)
can be used to obtain an optimal bin size of a regular histogram ingeneral. Hence, by using Eq.C.3 to estimate the sampling error termin Eq.3.13, we obtain
(
)
N
N
X
X
2
1
1
1
1
CnF IX (∆) =ki −k2n∆2 n − 1 N i=1n (n − 1) N i=1 iµ X ¶2
1 X 2
1
−
ki −
ki
N
N
N
X
2
n+1
1
=
− 2
(C.5)
ki2 − 2 .
(n − 1) ∆Tn (n − 1) ∆T i=1
TThis cost function is identified with Rudemo’s estimator for his riskfunction (Rudemo, 1982).
The relations between the cost function of the samples derived fromthe Poisson point process and CnF IX (∆) is given by
Cn (∆) =1 X 2
1n − 1 F IX
Cn (∆) + 3
ki +
.
n
n ∆T
nT 2

(C.6)

Due to the second term, the expected optimal bin size of a Poisson
point process is slightly larger than that obtained under fixed sample
sizes. This result arises from the fact that the histograms constructedfrom the data derived from a Poisson point process contains an additional fluctuation in its total sample size compared to the histograms
constructed from the data whose size is fixed. However, since the difference is on the order of O (n−3 ), the two cost functions yield indistinguishable optimal bin size in practice.

50

C.2

A line-graph histogram

From Eq.B.9, the cost function for a line-graph histogram becomes½¾
1 (+, +) 1 (−, −) 1 (+, −)
1
c
+ c
+ c
Cn (∆) =3
3(n∆)2 3©¡¢¡¢ª1(+, 0)(−, 0)(+, ∗)(−, ∗)
−
c
+
c
+
c
+
c
(n∆)2
¢ ¡ (+, ∗)¢ª1 ©¡ (+, 0)(−, 0)(−, ∗)
+
c̄
+
c̄
−
c̄
+
c̄
.
(C.7)
n∆2
The covariance with respect to the bins can be easily calculated asc(p,q) ≡N´³´
1 X ³ (p)(q)
ki − k̄ (p) ki − k̄ (q) ,N i=1

where p, q = {+, −, 0, ∗} and k̄ (p) = N1
cost function is obtained by usingc̄(+, 0) + c̄(−, 0) =PN(p)i=1 ki .

(C.8)

The last term in theNNo
1 1 X (0)1
1 X n (+)(−)
(0)
ki + kiki ,ki −n − 1 N i=1n (n − 1) N i=1(C.9)
c̄(+, ∗) + c̄(−, ∗) =NNo
1 1 X (∗)
1
1 X n (+)(−)(∗)
ki + kiki .
ki −n − 1 N i=1n (n − 1) N i=1(C.10)
These equations are derived from the argument below under theconstraint on the total sample size. We denote the covariance with(p, q)
respect to the sampled data in theith bin by c̄i , which is defined as
Ã
!Ã
!
n
(p)
(q)
1 X
ki
ki
(p, q)
(p)
(q)
c̄i
≡
ki (j) −ki (j) −.
(C.11)
n − 1 j=1nn51

The covariance averaged over the bins is given by c̄(p, q) = N1
(p, q)
where c̄i
is given by

PN(p, q),
i=1 c̄in´³´
1 X ³ (p)
(p)
(q)
(q)
ui (j) − ūiui (j) − ūin − 1 j=1

(p,q)
c̄i
=

n
n
1 X (p)(q)(p) (q)ui (j)ui (j)−ū ū .
n − 1 j=1
n−1 i i=(C.12)

Since the jth sample that enters the ith bin locates in either a left side[−∆/2, 0] or a right side [0, ∆/2] of the bin, we have
n
X

(+)
(0)
ui (j)ui (j) +

n
X

j=1

j=1

n
X

n
X

(+)

(∗)

ui (j)ui (j) +

j=1

(−)

(0)

(0)

(C.13)

(−)

(∗)

(∗)

(C.14)ui (j)ui (j)=ki ,ui (j)ui (j)=ki .

j=1
(p)

These equations hold because ui (j) = {0, 1}. Note that(p) (q)(p) (q)ūi ūi =ki ki /n2 .

(C.15)

By using these equations we obtainno11(+) (0)(−) (0)(0)ki ki + ki ki
,
ki −
n−1
n (n − 1)no11(+, ∗)(−, ∗)(∗)(+) (∗)(−) (∗)
c̄i
+ c̄i
=
k −
k k + ki ki.
n−1 in (n − 1) i i(+, 0)

c̄i(−, 0)

+ c̄i

=

52

(C.16)(C.17)

References
Blahut, R. E. (1987). Principles and practice of information theory.
Reading, Mass: Addison-Wesley.
Cover, T. M., & Thomas, J. A. (1991). Elements of information theory.
New York: John Wiley & Sons, Inc.
Cox, R. D. (1962). Renewal theory. London: John Wiley & Sons Inc.
Daley, D., & Vere-Jones, D. (1988). An introduction to the theory ofpoint processes. New York, USA: Springer-Verlag.
Devroye, L., & Györfi, L. (1985). Nonparametric density estimation :the l[1] view. New York: Wiley.
Diggle, P. J. (1983). Statistical analysis of spatial point patterns. London: Academic Press Inc.
Freedman, D., & Diaconis, P. (1981). On the histogram as a density
estimator - l2 theory. Zeitschrift Fur Wahrscheinlichkeitstheorie
Und Verwandte Gebiete, 57 (4), 453-476.
Gardiner, C. (1985). Handbook of stochastic methods: for physics,
chemistry and the naural sciences. Berlin: Springer-Verlag.
Hall, P. (1990). Akaikes information criterion and kullback-leibler loss
for histogram density-estimation. Probability Theory and Related
Fields, 85 (4), 449-467.
Izenman, A. J. (1991). Recent developments in nonparametric densityestimation. Journal of the American Statistical Association,
86 (413), 205-224.
Kanazawa, Y. (1993). Hellinger distance and akaike information criterion for the histogram. Statistics & Probability Letters, 17 (4),293-298.
53

Koyama, S., & Shinomoto, S. (2004). Histogram bin width selection for time-dependent poisson processes. Journal of Physics
A-Mathematical and General, 37 (29), 7255-7265.
Kubo, M., R. Toda, & Hashitsume, N. (1985). Nonequilibrium statistical mechanics. Berlin: Springer-Verlag.
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.
Scott, D. W. (1979). Optimal and data-based histograms. Biometrika,
66 (3), 605-610.
Shimazaki, H., & Shinomoto, S. (2007a). A method for selecting thebin size of a time histogram. Neural Computation. (in Press)
Shimazaki, H., & Shinomoto, S. (2007b). A recipe for optimizing
a time-histogram. In B. Schölkopf, J. C. Platt, & T. Hofmann(Eds.), Advances in neural information processing systems 19.
Cambridge, MA: MIT Press. (in Press)
Snyder, D. (1975). Random point processes. New York: John Wiley &Sons, Inc.
Taylor, C. C. (1987). Akaike information criterion and the histogram.
Biometrika, 74 (3), 636-639.
Wand, M. P. (1997). Data-based choice of histogram bin width. American Statistician, 51 (1), 59-64.

54