Biometrika (1979), 66, 3, pp.

605-10 605
Printed in Great Britain

On optimal and data-based histograms

BY DAVID W. SCOTT
Department of Mathematical Sciences, Rice University, Houston, Texas

SUMMARY
In this paper the formula for the optimal histogram bin width is derived which asymptotic-
ally minimizes the integrated mean squared error. Monte Carlo methods are used to verify
the usefulness of this formula for small samples. A data-based procedure for choosing the bin
width parameter is proposed, which assumes a Gaussian reference standard and requires only
the sample size and an estimate of the standard deviation. The sensitivity of the procedure is
investigated using several probability models which violate the Gaussian assumption.

Downloaded from http://biomet.oxfordjournals.org/ at Russian Archive on December 26, 2013

Some key words: Frequency distribution; Histogram; Nonparametrio density estimation; Optimal bin
width.

1. INTRODUCTION
The histogram is the classical nonparametric density estimator, probably dating from the
mortality studies of John Graunt in 1662 (Westergaard, 1968, p. 22). Today the histogram
remains an important statistical tool for displaying and summarizing data. In addition it
provides a consistent estimate of the true underlying probability density function. Present
guidelines for constructing histograms do not directly address the issues of estimation bias
and variance. Rather, they draw heavily on the investigator's intuition and past experience.
In this paper we propose new guidelines that reduce the subjectivity involved in histogram
construction by considering a mean squared error criterion.

2. BACKGROUND
We consider only histograms denned on an equally spaced mesh {tni; — oo < i < oo} with bin
width hn = tnH+1) — tni, where n denotes the sample size and emphasizes the dependence of
the mesh and bin width on the sample size. For a fixed point x, the mean squared error of a
histogram estimate,/(x), of the true density value,/(a;), is defined by
MSE(Z) =
For a random sample of size n from /, Cencov (1962) proved that MSE(X) asymptotically
converges to zero at a rate proportional to n~%ls, that is, MSE (a;) = (^(n"278). This rate is fairly
close to the Cram6r-Rao lower bound of 0{n~x). The integrated mean squared error repre-
sents a global error measure of a histogram estimate and is defined by

IMSE= [E{f{x)-f{x)Ydx.

Since it is the shape of the density that is of most interest, the IMSE is more relevant than the
mean squared error of the density height. The IMSE of a histogram also converges to zero
as 0(n- 2/s ).
To achieve these rates of convergence requires proper choice of the two parameters of the
histogram, the bin width hn and the relative position of the mesh. The latter is determined by
 606 D A V I D W. SCOTT

any particular mesh point, say tn0. Statistical texts suggest various methods for choosing
these two parameters. First the bin width is determined indirectly by choosing an appropriate
number of bins over the sample range. Most authors advise that 6-20 bins are usually
adequate for real data sets (Haber & Runyon, 1969, p. 33; Guttman & WUks, 1965, p. 59).
Larson (1975, p. 15) suggests using 1 + 2-2log10n bins as a first choice, similar to a formula
proposed by Sturges in 1926. The final choice for hn is a convenient whole number or fraction,
often related to the accuracy with which the data are measured. Next, tn0 is picked so that
the data do not fall on the bin boundaries. If we assume that the data are measured to
infinite accuracy, then the choice of tn0 becomes less important as the sample size increases.
Since we are focusing on consistency, we shall assume tnQ = 0 in the sequel. However, the
choice of hn is quite important. If hn is too small, then the histogram will be too rough; on
the other hand, if hn is too large, then the histogram will be too smooth, equivalent statistic-
ally to large variance and large bias, respectively. The proper choice for hn should balance
the bias and variance by minimizing, for example, the integrated mean squared error.

Downloaded from http://biomet.oxfordjournals.org/ at Russian Archive on December 26, 2013

In the past 20 years new nonparametric density estimators have been proposed and
investigated (Tapia & Thompson, 1978; Wegman, 1972). The most extensively treated of
these new estimators is the kernel probability density estimator developed by Rosenblatt
(1956) and Parzen (1962). The kernel estimator is also consistent but with IMSB = O(n~A/i),
an improvement over the histogram. In spite of these advances, the histogram will almost
surely retain important roles in data representation and density estimation, since it is simple
to compute and easily understood. Fortunately, by using techniques employed in kernel
density estimation consistency proofs, it is now possible to derive the optimal choice for the
bin width hn of a histogram.

3. DEBIVATION OF THB OPTIMAL HISTOGRAM BIN WIDTH

Suppose that a^,..., xn is a random sample from a continuous probability density function
/ with two continuous and bounded derivatives. We shall need to identify the bin interval
that contains a fixed point x as n varies. Let In(x) be that interval and let tn(x) denote the
left-hand endpoint of IJ^x). Define the bin probability

For y in In(x) we have, using Taylor's expansion,/(y) =f(x)+f'(x)(y—z) + O(A*). Therefore

C
J

Let vn(x) be the number of values falling in In(x). Then vn(x) has a binomial distribution
B{n,pn(x)}. The histogram estimate is given by the random variable
/(*) = vn(x)l{nhn),
with expectation

= /(*) + l* n /'(*) -/'(*) {x - tn(x)} + 0{h*).

Therefore the bias is
 On optimal and data-based histograms 607
Now the variance of the histogram estimate at a; is given by

f
Combining, we have that
MSE (x) = f{x)l(nhn) + \h\f\xf +/'(*)*{* - tn{x)f - hn /'(*)»{* - tn(x)} + O(l/n + A*). (1)
Integration of equation (1) over the real line implies that

IMSE = 1/(T*AJ + *A» J7'(*)» dx + jf'izWx - tn(x)}* dx . •

(as)«{*-U*)}*«' + O(l/n + A»). (2)

Downloaded from http://biomet.oxfordjournals.org/ at Russian Archive on December 26, 2013

Recall that {tni} denotes the mesh. Then the third term in equation (2) may be written as

i—cojtrt i—oo JO

by a change of variables. Now/'(< ni +y) =/'(^ B i) + ^(^n)> ^ ^ a * (3) becomes

JO J-co

by standard numerical integration approximations. A similar analysis for the fourth term
in (2) yields

Therefore

IMSE = ll(nhn)+&hl f00 /'(*)«eto + O(l/n + A«). (4)

J —00

Minimizing the first two terms in (4), we obtain

which, asymptotically, is the optimal choice for An.

We can estimate how the IMSE changes for poor choices of the bin width by using (4).
For any density and any positive constant c, the IMSE using the bin width ch* is larger than
the minimum IMSE by the factor (cs + 2)/(3c). Thus a bin width 50% too small implies an
IMSE 42% too large. We remark that a change of scale in the density function results in a
similar scaling of the optimal hn since y = x\a leads to J f'(y)ady = a 3 / f'(x)%dx. For Gaussian
data, h* = 2x 3 1/8 Tr1-"1 cm-"3.

4. SMAIii SAMPLE PBOPEBTIES

The formula for h* is based on an asymptotic expression. To investigate the small sample
properties of the IMSE, we undertook a fairly extensive Monte Carlo study of standardized
Gaussian data. For a range of values of A, the integrated squared error was computed exactly
608 DAVID W. SCOTT

for each of 1000 generated samples and then averaged over the number of repetitions to
obtain an estimate of the IMSE. The optimal bin widths predicted by equation (5) were quite
close to the empirically observed optimal bin widths for the Monte Carlo study even for
samples as small as 25. The estimated IMSE also increased as (c3 + 2)/(3c) for bin widths
differing from the empirically optimal bin width by the factor c.

5. DATA-BASED HISTOGRAMS
The optimal choice for hn requires knowledge of the true underlying density/. This know-
ledge is rare. In another context Tukey (1977, p. 623) has suggested using the Gaussian
density as a reference standard, to be used cautiously but frequently. Therefore, we propose
the data-based choice for the bin width
hn = 3-49an-1/8, (6)
where a is an estimate of the standard deviation. Although the Gaussian density forms the

Downloaded from http://biomet.oxfordjournals.org/ at Russian Archive on December 26, 2013

basis of (6), this assumption is not so strong as a parametric Gaussian assumption, i.e. use of
equation (6) on non-Gaussian data will not result in a histogram that looks Gaussian. For
density functions with equal variances, the data-based choice (6) results in the same bin
width. To show that (6) is useful for a large class of densities, we considered Gaussian and
non-Gaussian densities with equal variances and observed how their theoretically optimal
bin widths (5) differed. In particular we considered three models of non-Gaussian behaviour:
skewed, heavy-tailed and bimodal densities.
As a model of skewed data, we used a log normal density with variance equal to
«^(«; a -l)exp(-2y/8) and skewness (tu* + 2){wi- 1)*, where w = exp(^-*). In Fig. l(a) we
plot the ratio of the optimal bin width for the log normal density to the optimal bin width
for a Gaussian density with the same variance as a function of the log normal skewness.
This ratio does not depend on the sample size. We see that using the Gaussian h* over-
smoothes a log normal density; however, for skewnesses as great as one, the difference is
less than 30%. A similar plot results when using a gamma probability density model.
We used Student's tr density to model heavy-tailed data. The variance and kurtosis are
r/(r —2) and 6/(r —4), respectively. In Fig. l(b) we plot the ratio of bin widths as a function
of the kurtosis, connecting the discrete points by a solid line for convenience. The insensitivity
of the data-based choice for hn for any moderate kurtosis is apparent.
(a) Log normal density (b) Student's density (o) Mixture of two Geussianfl

0-5

6 6
Skewnees coefficient KurtosU coefficient Distance between modes
Fig. 1. Ratio of theoretical bin width for several non-Gaussian probability densities to the
theoretical bin width for a Gaussian density with the same variance.

Finally as a model of bimodal data, we used a mixture of Gaussian distributions,

$ ( / * , l) + pT(/n, 1), with variance l + /x2. In Fig. l(c) we give a similar ratio of theoretical
bin widths as a function of 2^,, the distance between the two modes. For strongly bimodal
data (ft, > 1-5), the ratio falls below 0-8, corresponding to oversmoothing of the bimodal data.
When distinctly bimodal data are encountered, the data-based histogram may be inadequate.
 On optimal and data-based histograms 609
Thus the data-based algorithm leads to An's that are generally too big for all our models of
non-Gaussian data. A correction factor may be applied to the data-based hn by computing
the sample skewness or kurtosis and reading the correction factor from Fig. 1. A histogram
does not exhibit sensitivity to small changes from the optimal hn as is evident from the
discussion after (5). We do not advocate using exactly the hn suggested by (6), but rather a
convenient choice either slightly larger or smaller.

6. EXAMPLES
In Fig. 2 we display three histograms of a Monte Carlo N(0,1) sample of size 1000 which
has a sample standard deviation equal to 1-011 with h = 0-176, 0-353 and 0-706, the second
choice obtained from (6). Many statisticians prefer a smaller bin width and a rougher histo-
gram, leaving the final smoothing to be done by eye.

Downloaded from http://biomet.oxfordjournals.org/ at Russian Archive on December 26, 2013

0-176 (b) h = 0-363 (o) h = 0-706

- 3 - 2 - 1 0 1 2 - 3 - 2 - 1 0 1 2 - 3 - 2 - 1 0 1 2
•Scale of observations Scale of observations Scale of observations

Fig. 2. Histograms of 1000 pseudorandom Gaussian numbers for three bin widths: the data-based
choice and that choice perturbed by a factor of 2.

To illustrate extremely large sample sizes, Kendall & Stuart (1969, p. 8) consider a histo-
gram of the ages of 301,785 Australian bridegrooms (1907-14) with a bin width of 3 years.
The sample standard deviation and skewness for these data are 7-97 and 1-93, respectively.
Thus the data-based choice for A is 0-41 years. Applying a skewness correction factor of
0-43 using Fig. l(a), the final data-based choice is 0-18 years. Thus the sample is of sufficient
size to use a bin width of 1 year or even 3 months if the data were recorded to sufficient
accuracy.

7. DISCUSSION
We have considered the optimal construction of histograms given either knowledge of the
true underlying density or, more commonly, given only the data. Waterman & Whiteman
(1978) have recently carried out a similar attack for Rosenblatt's kernel estimator. Kernel
estimates converge faster than histograms to the true density, and therefore integrated mean
squared error is more sensitive to the choice of the smoothing parameter; see also Silverman
(1978). Furthermore, kernel estimates require the entire data set for evaluation. Thus in
some modern automated data collectors, it is often more economical to summarize sequenti-
ally relatively more samples, calibrating the histogram using a small training sample.
Some recently developed nonparametric techniques for density estimation start with a
histogram and then smooth it; see, for example, Boneva, Kendall & Stefanov (1971). Our
procedures could be used to construct the required histogram directly from the data. We
remark that our analysis extends easily to histograms in higher dimensions.
34
610 DAVID W. SCOTT

It should be possible to further reduce the integrated mean squared error by using an
unequally spaced mesh. However, the algorithms required would surely be iterative and
would require the entire data set. It is easier to discount rougher estimates in the tails or to
construct a rootgram as suggested by Tukey (1977, p. 543).

This research was supported in part by the National Heart, Lung, and Blood Institute, the
National Institutes of Health, the Department of Health, Education and Welfare. The
author would like to thank a referee for helpful commente.

REFERENCES
BONEVA, L. I., KENDALL, D. G. & STEFANOV, I. (1971). Spline transformations: Three new diagnostic aids
for the statistical data-analyst (with discussion). J. R. Statist. Soc. B 83, 1-70.
CENOOV, N. N. (1962). Estimation of an unknown distribution density from observations. Soviet Math. 3,
1559-62.
GUTTMAN, I. & WrLKS, S. S. (1965). Introductory Engineering Statistics. New York: Wiley.

Downloaded from http://biomet.oxfordjournals.org/ at Russian Archive on December 26, 2013

HABKR, A. & RUNYON, R. P. (1969). General Statistics. Reading, MJUSB: Addison-Wesley.
KENDALL, M. G. & STUART, A. (1969). The Advanced Theory of Statistics, Vol. 1. London: Griffin.
LABSON, H. J. (1975). Statistics: An Introduction. New York: Wiley.
PABZEN, E. (1962). On estimation of a probability density function and mode. Ann. Math. Statist. 38,
1065-76.
ROSENBLATT, M. (1956). Remarks on some nonparamteric estimates of a density function. Ann. Math.
Statist. 27, 832-7.
SELVEBMAN, B. W. (1978). Choosing the window width when estimating a density. Biometrika 65, 1-11.
STUBGEB, H. A. (1926). The choice of a class interval. J. Am. Statist. Assoc. 21, 65-6.
TAPIA, R. A. <fc THOMPSON, J. R. (1978). Nonparametric Probability Density Estimation. Baltimore:
Johns Hopkins.
TUXBY, J. W. (1977). Exploratory Data Analysis. Reading, Mass: Addison-Wesley.
WATERMAN, M. S. <fc WHITEMAN, D. E. (1978). Estimation of probability densities by empirical density
functions. Int. J. Math. Educ. Sd. Technol. 9, 127-37.
WEGMAN, E. J. (1972). Nonparametric probability density estimation: I. A summary of available
methods. Technometrics 14, 533-46.
WESTERGAABD, H. (1968). Contributions to the History of Statistics. New York: Agathon.

[Received February 1979. Revised May 1979]