This page INSTANTLY generates an optmized histogram of YOUR DATA.

*Solutions may be found at FAQ.
*The data is processed on your own computer (the data is NOT transferred to our server).
*The program searches number of bins from 2 to 200.
*This program does not work in Opera. Mac users may experience an extremely slow response.

07/06/16 Ver. 1.0
06/11/08 Ver. 0.10 Improved data sheet output.

I. Divide the data range into bins of width . Count the number of events that enter the i'th bin.

and

IV. Repeat i-iii while changing . Find that minimizes .

An Optimized Histogram

Data: the duration for eruptions of the Old Faithful geyser
in Yellowstone National Park (in minutes)

From the observed data only, the method estimates a binwidth that minimizes expected L2 loss between the histogram and an unknown underlying density function. An assumption made here is merely that samples are drawn from the density independently each other.

Download a Full paper
Shimazaki H. and Shinomoto S., A method for selecting the bin size of a time histogram Neural Computation (2007) Vol. 19(6), 1503-1527

Short paper
Shimazaki H. and Shinomoto S., A recipe for optimizing a time-histogram Neural Information Processing Systems (2007) Vol. 19, 1289-1296

Presentation slides and FAQ

New original paper on Kernel Density Estimation
Shimazaki H. and Shinomoto S., Kernel Bandwidth Optimization in Spike Rate Estimation. Journal of Computational Neuroscience (2009) doi:10.1007/s10827-009-0180-4 [Open Access: PDF, LINK]

Go to Kernel Bandwidth Optimization page.

1. An unknown underlying rate that generates observable events.
2. Event sequences generated from the unknown rate.
3. Histograms constructed from the events. It is not obvious which bin size should be used.
4. Mean Integrated Squared Error (MISE), a measure of the goodness-of-the-fit of a histogram to the unknown rate. A histogram with the bin size that minimizes MISE is optimal.
5. Note that we can not directly compute the MISE since we do not know the underlying rate. However, the MISE can be estimated from the data.
6. The optimal bin size can be estimated as the one that minimizes the estimated MISE.

• If you choose too small bin size, the bar height at each bin suffers from significantly large statistical fluctuation due to the paucity of samples in each bin.
• If you choose too large bin size, the histogram can not represent the shape of the underlying distribution because the resolution isn't good enough.

There is an optimal bin size that remedies these two opposite problems. The Java Applet below demonstrates the problem of the bin size selection.

In the applet, data is shown in the upper box and its histogram in the bottom box. The data is an example from Neuroscience: a timing of neuronal firing (spikes) obtained by repeating trials. Twenty sequences (trials) are drawn in the data box. However, the data can be anything (weight, height, or test score etc.), and can be one sequence.

The histogram is shown in the bottom box by red color. The blue line is the distribution (or rate) that produce data samples in the box above. Our aim is to choose a histogram that best represents this blue line.

As a first step, change the bin size of a histogram in the applet by using the scroll bar at the bottom. With too small a bin size, you get a jagged histogram. With too large a bin size, you obtain a flat histogram.

Now bring the scroll bar left to make the bin size thin, and push `redraw' button several times. The `redraw' button regenerates sample points based on the blue density distribution. Please confirm that the shape of the histogram drastically changes. If you are told to tell which of the two hills of the distribution is taller, can you tell? Although the right hill is taller than the left (blue line), the maximum height of a histogram appears at left side quite often.

As a second step, bring your scroll bar right to make the bin size broad, and push `redraw' several times. Please confirm that the shape of the histogram does not change very much. In addition, you may notice the right hand side of the hills is higher than the left hand side most of the time. However, can you tell where the highest point of the hill is? The resolution of the histogram is not good enough.

When you make a histogram, you need to choose a bin size that compromises the conflict between sampling error and resolution. Now please check the `error check box' on. Yellow area appears, which indicates an error between the histogram and the underlying distribution. You may turn off the histogram check box. Please examine where a total error becomes the smallest by your eyes. You would find that the error becomes the smallest when you bring the scroll bar handle about one fourth of the total length from the left. Please push the `redraw' button several times to check the statistical fluctuation of errors. The optimal bin size of a histogram is the bin size that produces the smallest total errors (yellow area) on average of many realizations of sampled data (i.e. average over many pushes of `redraw' button).

At this moment, you might say, "Okay, I now understand what an optimal bin size is. But, how can we know it? After all, we do not know the blue underlying distribution. So we can not know the errors. We can never know the optimal bin size from data."

This statement is not true. Indeed, the errors can be estimated from the data. The above method computes the estimated error for several candidate bin sizes, and choose the bin size that produces the smallest `estimated' error. For the details of the method, please refer to Shimazaki and Shinomoto in Neural Computation.

A. Yes.

Q. I obtained a very small bin width, which is likely to be erroneous. Why?

A. You probably searched bin widths smaller than the sampling resolution of your data. Please begin the search from the bin width as large as your sampling resolution. In addition, please check if your data contains duplicated samples. A small bin width will be selected for such highly correlated samples. See below for the solution to the duplicated data due to the low sampling resolution.

Q. I have data with low sampling resolution.

A. To obtain the optimal bin size correctly, it is recommended to replace your data x with x+r, where r is an uniform random variable drawn from [-dx/2, dx/2]. dx is sampling resolution of data acquisition. You repeat this to all samples x with independently drawn r. Use this randomized data to obtain the optimal bin size. To draw histogram, use the original data. See the difference with the following examples.

Copy and paste rounded data

Q. Can I apply the method to the data composed of integer values?

A. Yes. Please note that the resolution of your data is an integer. It is recommended to search only the bin widths which are a multiple of the resolution of your sampled data. Do not search the bin width smaller than the resolution.

(The current version of web application can NOT be used for computing only integer bin widths.)

Q. I want to make a 2-dimensional histogram. Can I use this method?

A. Yes. The 'bin size' of a 2-d histogram is the area of a segmented square cell. The mean and the variance are simply computed from the event counts in all the bins of the 2-dimensional histogram.

Matlab demo program for selecting bin size of 2-d histogram.

(The current version of web application can NOT be used for computing 2-dimensional histogram.)

See also 2-d kernel density estimation.

Q. Can I use unbiased variance as the variance in the formula of bin size selection?

A. No. Please use the (biased) variance displayed in the method. Please note that a built-in function for variance calculation in a software may returns unbiased variance in default.

A. They should be different. Three assumptions were made to obtain the Scott's result. First, the Scott's result is asymptotically true (i.e. it is true for large sample size n). Second, the scaling exponent -1/3 is true if the density is a smooth function. Third, the coefficient 3.49 was obtained, assuming the Gauss density function as a reference. The present method does not require these assumptions.

Q.What assumption was made in your method?

A. The method only assumes that the data are sampled independently each other (an assumption of a Poisson point process). No assumption was made for the underlying density function (for example, unimodality, continuity, existence of derivatives, etc.).

Q. Is the assumption in this method distinct from that in other methods?

A. Yes. In classical density estimation, the sampling size is fixed. Here, the total data size is not fixed, but assumed to obey a Poisson distribution. Note that the MISE is ensemble performance of the histograms with 'many' realizations of the experiment. The histograms constructed from a Poisson point events have more variability than the histograms constructed from samples with a fixed amount. This leads to a choice of wider optimal bin size for a histogram under a Poisson point process assumption. However, this difference becomes negligible as we increase the data size.

Q. You provide the method for optimizing kernel density estimate, too. Which of the methods, a histogram or kernel, do you recommend for density estimation?

A. We generally recommend to use the kernel density estimate. Please check the kernel density optimization page, too.

A bar-graph histogram is constructed simply by counting the number of events that belong to each bin of width . We divide the data range into intervals with each length given by . The number of events accumulated from all sequences in the th interval is counted as . The bar height at the th bin is given as .

In this thesis, we assess the goodness of the fit of the estimator to the underlying rate over the data range by the mean integrated squared error (MISE),

MISE

(3.3)

By segmenting the range into intervals of size , the MISE defined in Eq.(3.3) can be rewritten as

MISE

(3.4)

MISE

(3.5)

MISE
		(3.6)

The second term of Eq.(3.5) can further be decomposed into two parts,

(3.7)

(3.8)

	MISE
		(3.9)

		(3.10)

(3.11)

(3.14)

(3.15)

(3.17)

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

The matlab function sshist.m is now available. [2008/11/19]
To plot an optimized histogram of the data x (a vector contains samples as an element), type

optN = sshist(x); hist(x,optN);

Below is a simple sample program sample.m .

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2006 Author Hideaki Shimazaki
% Department of Physics, Kyoto University
% shimazaki at ton.scphys.kyoto-u.ac.jp
% Please feel free to use/modify this program.
%
% Data: the duration for eruptions of
% the Old Faithful geyser in Yellowstone National Park (in minutes)
clear all;
x = [4.37 3.87 4.00 4.03 3.50 4.08 2.25 4.70 1.73 4.93 1.73 4.62 ...
     3.43 4.25 1.68 3.92 3.68 3.10 4.03 1.77 4.08 1.75 3.20 1.85 ...
     4.62 1.97 4.50 3.92 4.35 2.33 3.83 1.88 4.60 1.80 4.73 1.77 ...
     4.57 1.85 3.52 4.00 3.70 3.72 4.25 3.58 3.80 3.77 3.75 2.50 ...
     4.50 4.10 3.70 3.80 3.43 4.00 2.27 4.40 4.05 4.25 3.33 2.00 ...
     4.33 2.93 4.58 1.90 3.58 3.73 3.73 1.82 4.63 3.50 4.00 3.67 ...
     1.67 4.60 1.67 4.00 1.80 4.42 1.90 4.63 2.93 3.50 1.97 4.28 ...
     1.83 4.13 1.83 4.65 4.20 3.93 4.33 1.83 4.53 2.03 4.18 4.43 ...
     4.07 4.13 3.95 4.10 2.27 4.58 1.90 4.50 1.95 4.83 4.12];

x_min = min(x);
x_max = max(x);

N_MIN = 4;              % Minimum number of bins (integer)
                        % N_MIN must be more than 1 (N_MIN > 1).
N_MAX = 50;             % Maximum number of bins (integer)

N = N_MIN:N_MAX;                      % # of Bins
D = (x_max - x_min) ./ N;             % Bin Size Vector


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Computation of the Cost Function
for i = 1: length(N)
	edges = linspace(x_min,x_max,N(i)+1);	% Bin edges

	ki = histc(x,edges);            % Count # of events in bins
	ki = ki(1:end-1);

	k = mean(ki);                   % Mean of event count
	v = sum( (ki-k).^2 )/N(i);      % Variance of event count

	C(i) = ( 2*k - v ) / D(i)^2;    % The Cost Function

end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Optimal Bin Size Selectioin
[Cmin idx] = min(C);
optD = D(idx);                         % *Optimal bin size
edges = linspace(x_min,x_max,N(idx)+1);  % Optimal segmentation


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Display an Optimal Histogram and the Cost Function
subplot(1,2,1); hist(x,edges); axis square;
subplot(1,2,2); plot(D,C,'k.',optD,Cmin,'r*'); axis square;

Copy and Paste a sample program sample.R, then run sshist(faithful[,1]) .

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# 2006 Author Hideaki Shimazaki
# Department of Physics, Kyoto University
# shimazaki at ton.scphys.kyoto-u.ac.jp
# Please feel free to use/modify this program.
sshist <- function(x){

	N <- 2: 100
	C <- numeric(length(N))
	D <- C

	for (i in 1:length(N)) {
		D[i] <- diff(range(x))/N[i]

		edges = seq(min(x),max(x),length=N[i])
		hp <- hist(x, breaks = edges, plot=FALSE )
		ki <- hp$counts

		k <- mean(ki)
		v <- sum((ki-k)^2)/N[i]

		C[i] <- (2*k-v)/D[i]^2	#Cost Function
	}

	idx <- which.min(C)
	optD <- D[idx]

	edges <- seq(min(x),max(x),length=N[idx])
	h = hist(x, breaks = edges )
	rug(x)

	return(h)
}

Copy and Paste a sample program sample.nb, then ctrl+shift.

\!\(<< \ Graphics`Graphics`\n
  << Statistics`DescriptiveStatistics`\[IndentingNewLine]
  \(x = {4.37, 3.87, 4.00\ , 4.03, \ 3.50, \ 4.08, \ 2.25, \ 4.70, \ 1.73, \ 
        4.93, \ 1.73, \ 4.62, \ 3.43, \ 4.25, \ 1.68, \ 3.92, \ 3.68, \ 
        3.10, \ 4.03, \ 1.77, \ 4.08, \ 1.75, \ 3.20, \ 1.85, 4.62, \ 1.97, \ 
        4.50, \ 3.92, \ 4.35, \ 2.33, \ 3.83\ , 1.88, \ 4.60, \ 1.80, \ 
        4.73, \ 1.77, \ 4.57, \ 1.85, \ 3.52, \ 4.00, \ 3.70, \ 3.72, \ 
        4.25, \ 3.58, \ 3.80, \ 3.77, \ 3.75, \ 2.50, 4.50, \ 4.10, \ 3.70, \ 
        3.80, \ 3.43, \ 4.00, \ 2.27, \ 4.40, \ 4.05, \ 4.25, \ 3.33, \ 
        2.00, \ 4.33, \ 2.93, \ 4.58, \ 1.90, \ 3.58, \ 3.73, \ 3.73, \ 
        1.82, \ 4.63, \ 3.50, \ 4.00, \ 3.67, \ 1.67, \ 4.60, \ 1.67, \ 
        4.00, \ 1.80, \ 4.42, \ 1.90\ , 4.63, \ 2.93, \ 3.50, \ 1.97, \ 
        4.28, \ 1.83, \ 4.13, \ 1.83, \ 4.65, \ 4.20, \ 3.93, \ 4.33, \ 
        1.83, \ 4.53, \ 2.03, \ 4.18, \ 4.43, \ 4.07, \ 4.13, \ 3.95, \ 
        4.10, \ 2.27, \ 4.58, \ 1.90, \ 4.50, \ 1.95, \ 4.83, \ 
        4.12};\)\[IndentingNewLine]
  c = {}; Δ = {}; \ d = {}; p = {};\[IndentingNewLine]
  For[i = 1, 
    i < 50, \(i++\), \[IndentingNewLine]n = i + 3; \[IndentingNewLine]d = 
      Append[d, \((Max[x] - Min[x])\)/n]; \[IndentingNewLine]ki = 
      BinCounts[
        x, {Min[x], Max[x], 
          d[\([i]\)]}]; \[IndentingNewLine] (*Cost\ Function\
*) \[IndentingNewLine]k\  = \ Mean[ki]; \[IndentingNewLine]v\  = \ 
      VarianceMLE[ki]; \[IndentingNewLine]c = 
      Append[c, \((2  k - v)\)/d[\([i]\)]\^2]; \[IndentingNewLine]p = 
      Append[p, {d[\([i]\)], 
          c[\([i]\)]}];\[IndentingNewLine]]\[IndentingNewLine] (*Optimal\ Bin\
\ Size\ Selection*) \[IndentingNewLine]
  \(idx = \(Ordering[c]\)[\([1]\)];\)\[IndentingNewLine]
  \(d\^*\) = d[\([idx]\)]\[IndentingNewLine]
  \(edges = {};\)\[IndentingNewLine]
  \(edges = 
      Do[Append[edges, Min[x] + \(d\^*\)*i], {i, 0, 
          idx + 3}];\)\[IndentingNewLine] (*Display\ an\ Optimized\ Histogram\
*) \[IndentingNewLine]
  \(Histogram[x, \ HistogramRange -> {Min[x], Max[x]}, 
      HistogramCategories -> edges];\)\[IndentingNewLine]
  \(ListPlot[p];\)\)

Other applications for analyzing spike data: SULAB ( Prof. Shinomoto )