Optimal bin width for two dimensional histogram

Question

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?

Answer 1

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.

Answer 2

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.