# Optimal bin width for two dimensional histogram

There are lots of rules for selecting an optimal bin width in a 1D histogram (see for example)

I'm looking for a rule that applies the selection of optimal equal-bin widths on two-dimensional histograms.

Is there such a rule? Perhaps one of the well-known rules for 1D histograms can be easily adapted, if so, could you give some minimal details on how to do so?

• Optimal for what purpose? Please also note that 2D histograms will suffer from the same problems seen in ordinary histograms so you might want to turn your attention to alternatives such as kernel density estimates.
– whuber
Sep 5, 2014 at 22:47
• Is there a reason why you would not adapt something something simple like the $\sqrt(n)$ rule or Sturges' formula to your problem directly? Along each dimension you have the same number of readings anyway. If you want something a bit more sophisticated (eg. Freedman-Diaconis rule) you could "naively" take the max between the number of bins return for each dimension independently. You are essentially looking into a discretized (2d) KDE anyway so maybe that is your best choice anyway. Sep 5, 2014 at 22:51
• For the purpose of not having to chose a bin width manually hence subjectively? For selecting a width that will describe the underlying data with not too much noise and not too smoothed? I'm not sure I understand your question. Is "optimal" a too vague word? What other interpretations can you see here? How else could have I phrased the question? Yes, I am aware of KDE but I need a 2D histogram. Sep 5, 2014 at 22:51
• @usεr11852 Could you expand your comment in an answer, perhaps with some more details? Sep 5, 2014 at 22:53
• @Glen_b could you put that in the form of an answer? My knowledge of statistics is pretty limited and many of the things you say go right over my head, so as many details as possible would be appreciated. Sep 6, 2014 at 0:23

My advice would generally be that it's even more critical than in 1-D to smooth where possible i.e. to do something like kernel density estimation (or some other such method, like log-spline estimation), which tends to be substantially more efficient than using histograms. As whuber points out, it's quite possible to be fooled by the appearance of a histogram, especially with few bins and small to moderate sample sizes.

If you're trying to optimize mean integrated squared error (MISE), say, there are rules that apply in higher dimensions (the number of bins depends on the number of observations, the variance, the dimension, and the "shape"), for both kernel density estimation and histograms.

[Indeed many of the issues for one are also issues for the other, so some of the information in this wikipedia article will be relevant.]

This dependence on shape seems to imply that to choose optimally, you already need to know what you're plotting. However, if you're prepared to make some reasonable assumptions, you can use those (so for example, some people might say "approximately Gaussian"), or alternatively, you can use some form of "plug-in" estimator of the appropriate functional.

Wand, 1997$^{[1]}$ covers the 1-D case. If you're able to get that article, take a look as much of what's there is also relevant to the situation in higher dimensions (in so far as the kinds of analysis that are done). (It exists in working paper form on the internet if you don't have access to the journal.)

Analysis in higher dimensions is somewhat more complicated (in pretty much the same way it proceeds from 1-D to r-dimensions for kernel density estimation), but there's a term in the dimension that comes into the power of n.

Sec 3.4 Eqn 3.61 (p83) of Scott, 1992$^{[2]}$ gives the asymptotically optimal binwidth:

$h^∗=R(f_k)^{-1/2}\,\left(6\prod_{i=1}^dR(f_i)^{1/2}\right)^{1/(2+d)} n^{−1/(2+d)}$

where $R(f)=\int_{\mathfrak{R}^d} f(x)^2 dx$ is a roughness term (not the only one possible), and I believe $f_i$ is the derivative of $f$ with respect to the $i^\text{th}$ term in $x$.

So for 2D that suggests binwidths that shrink as $n^{−1/4}$.

In the case of independent normal variables, the approximate rule is $h_k^*\approx 3.5\sigma_k n^{−1/(2+d)}$, where $h_k$ is the binwidth in dimension $k$, the $*$ indicates the asymptotically optimal value, and $\sigma_k$ is the population standard deviation in dimension $k$.

For bivariate normal with correlation $\rho$, the binwidth is

$h_i^* = 3.504 \sigma_i(1-\rho^2)^{3/8}n^{-1/4}$

When the distribution is skewed, or heavy tailed, or multimodal, generally much smaller binwidths result; consequently the normal results would often be at best upper bounds on bindwith.

Of course, it's entirely possible you're not interested in mean integrated squared error, but in some other criterion.

[1]: Wand, M.P. (1997),
"Data-based choice of histogram bin width",
American Statistician 51, 59-64

[2]: Scott, D.W. (1992),
Multivariate Density Estimation: Theory, Practice, and Visualization,
John Wiley & Sons, Inc., Hoboken, NJ, USA.

• Glen, what is the $\sigma_k$ parameter in your answer? Sep 7, 2014 at 23:58
• My apologies; I've added the meaning as I understand it - population s.d. in dimension $k$. (I had a link to a pdf in a comment before, but after later figuring out it was a chapter from a book -- the Scott reference in my answer -- I removed the link.) Sep 8, 2014 at 0:29

Given you have a fixed number $N$ of data (ie. you have equal number of reading on both dimensions) you could immediately use:

1. The square-root rule rounded-down ($\sqrt{N}$), (ie. the Excel-way :) )
2. Sturges' rule ($\log_2N +1$),
3. Some other rule that is based only on the number of available data-points (eg. Rick's rule).

To find the common number of bins $M$ across each dimension.

On the hand, you might want to try something more robust like the Freedman–Diaconis rule which essentially defines the bandwidth $h$ as equal to:

$h= 2 IQR(x) N^{-1/3}$,

where IQR is the interquartile range of your data $x$. You then calculate the number of bins $M$ along each dimension as being equal to:

$M = \lceil (max(x)- min(x))/h \rceil$.

You do this across both dimensions of your data $x$; this gives you two, possibly different, numbers of bins that "should" be used across each dimension. You naively take the larger one so you do not "lose" information.

Yet, a fourth option would be to try to treat your sample as natively two-dimensional, calculate the norm for each of the sample points and then perform the Freedman–Diaconis rule on the sample's norms. ie.:

$x_{new} = \sqrt{x_1^2 + x_2^2}$

OK, here is some code and a plot for the procedures I describe:

rng(123,'twister');     % Fix random seed for reproducibility
N = 250;                % Number of points in our sample

A = random('normal',0,1,[N,2]);  % Generate a N-by-2 matrix with N(0,1)
A(:,2) = A(:,2) * 5;             % Make the second dimension more variable

% The sqrt(N) rule:
nbins_sqrtN = floor(sqrt(N));

% The Sturges formula:
nbins_str = ceil(log2(N) +1);

% The Freedman–Diaconis-like choice:
IQRs = iqr(A);              % Get the IQ ranges across each dimension
Hs = 2* IQRs* N^(-1/3);     % Get the bandwidths across each dimension
Ranges = range(A);          % Get the range of values across each dimension
% Get the suggested number of bins along each dimension
nbins_dim1 = ceil(Ranges(1)/Hs(1)); % 12 here
nbins_dim2 = ceil(Ranges(2)/Hs(2)); % 15 here
% Get the maximum of the two
nbins_fd_1 = max( [nbins_dim1, nbins_dim2]);

% The Freedman–Diaconis choice on the norms

Norms = sqrt(sum(A.^2,2));        % Get the norm of each point in th 2-D sample
H_norms = 2* iqr(Norms)* N^(-1/3);% Get the "norm" bandwidth
nbins_fd_2 = ceil(range(Norms)/ H_norms);   % Get number of bins

[nbins_sqrtN nbins_str nbins_fd_1 nbins_fd_2]

% Plot the results / Make bivariate histograms
% I use the hist3 function from MATLAB
figure(1);
subplot(2,2,1);
hist3(A,[ nbins_sqrtN nbins_sqrtN] );
title('Square Root rule');

subplot(2,2,2);
hist3(A,[ nbins_str nbins_str] );
title('Sturges formula rule');

subplot(2,2,3);
hist3(A,[ nbins_fd_1 nbins_fd_1]);
title('Freedman–Diaconis-like rule');

subplot(2,2,4);
hist3(A,[ nbins_fd_2 nbins_fd_2]);
title('Freedman–Diaconis rule on the norms');


As others have noted smoothing is almost certainly more appropriate for this case (ie. getting a KDE). I hope thought this gives you an idea about what I described in my comment regarding the direct generalization (with all the problems it might entail) of 1-D sample rules to 2-D sample rules. Notably, most procedure do assume some degree of "normality" in the sample. If you have a sample that clearly is not normally distributed (eg. it is leptokurtotic) these procedure (even in 1-D) would fail quite badly.

• I am glad you found it helpful Gabriel but I would propose accepting Glen_b 's answer. While convenient the methods described in my answer are heuristics generalizing from the 1-D case; I added them to illustrate the points I made in my comment. I truly appreciate the thumbs up (so you can still up-vote me as giving a useful answer if you haven't already :) ) but the correct answer is Glen_b's. Sep 6, 2014 at 20:58