2D Histogram Calculation

Quick computation of two-dimensional histogram of bivariate data

Description

function histmat = hist2(x, y, xedges, yedges)
 
  Extract 2D histogram data containing the number of events of [x , y] pairs that fall in each bin of the grid defined by xedges and yedges. The edges are vectors with monotonically non-decreasing values.

The code is optimized no loop inside, it can be useful in the case of large dataset.
 
 EXAMPLE
 
  events = 1000000;
  x1 = sqrt(0.05)*randn(events,1)-0.5; x2 = sqrt(0.05)*randn(events,1)+0.5;
  y1 = sqrt(0.05)*randn(events,1)+0.5; y2 = sqrt(0.05)*randn(events,1)-0.5;
  x= [x1;x2]; y = [y1;y2];
 
 For linearly spaced edges:
  xedges = linspace(-1,1,64); yedges = linspace(-1,1,64);
  histmat = hist2(x, y, xedges, yedges);
  figure; pcolor(xedges,yedges,histmat'); colorbar ; axis square tight;
 
 For nonlinearly spaced edges:
  xedges_ = logspace(0,log10(3),64)-2; yedges_ = linspace(-1,1,64);
  histmat_ = hist2(x, y, xedges_, yedges_);
  figure; pcolor(xedges_,yedges_,histmat_'); colorbar ; axis square tight;

• Acknowledgements: This submission has inspired the following:; hist2 for the people, histograms for ND data
• MATLAB release: MATLAB 6.5.1 (R13SP1)

[Nedialko Krouchev]

It's odd in itself that Matlab or stats don't seem to provide this function.
Laszlo has made good use of histc toward an even more useful purpose.

[Nedialko Krouchev]

The fix to the 'out of range values' is not good as it confounds the two tails of the distribution. Fixed it in my copy.

[Nedialko Krouchev]

The 'end-spiel'(main number processing, after dealing with zero bin indexes) part is overcomplicated and can be replaced by an one liner. Hardly as fast as claimed.
Please see the updated entry.

[Narupon Chattrapiban]

This code works fine for me
I like the way "histc" is use. Ingenious!

The out of range values should be dumped into a bin other than the bin that contains zeros.

The code might be rewritten as follow:
.
.
.

[xn, xbin] = histc(x,xedges);
[yn, ybin] = histc(y,yedges);

clear xn yn;

%xbin, ybin zero for out of range values
% (see the help of histc) force this event to
% inf
xbin(find(xbin == 0)) = inf; % -> [changed]
ybin(find(ybin == 0)) = inf; % -> [changed]

xnbin = length(xedges);
ynbin = length(yedges);

if xnbin >= ynbin
    xy = ybin*(xnbin) + xbin;
      indexshift = xnbin;
else
    xy = xbin*(ynbin) + ybin;
      indexshift = ynbin;
end

xyuni = unique(xy);
xyuni(end) = []; % -> [added]
hstres = histc(xy,xyuni);
clear xy;

histmat = zeros(ynbin,xnbin); % -> [changed]
histmat(xyuni-indexshift) = hstres;

[Nathan]

 

[Radford Juang]

Pretty useful. There's a slight bug though with how the final histmat is created.
1. It is transposed when xnbin >= ynbin.
2. It is incorrectly indexed when ynbin > xnbin (because there is a difference in the index being row major and column major)

[Zbigniew]

Matlab's hist3 calculates 2Dhistograms.

[John]

A surface can be created using

surf(xedges,yedges,histmat');

[John]

I'm having trouble with this routine. In the example below, the x and y ranges dont' seem to be correct, and also the distribution doesn't look right either.

function hist3_test
close all; clc;
events = 1000000;
xrange=0:.1:1;
nxs=length(xrange);
x=[]; y=[];
for xx=0:.1:1
    xs=xx+sqrt(0.005)*randn(round(events/nxs),1);
    if (xx<.3)
        ys=cos(xs)+ sqrt(0.001)*randn(round(events/nxs),1);
    else
        ys=cos(xs)+ sqrt(0.1)*randn(round(events/nxs),1);
    end
    x = [x; xs];
    y = [y; ys];
end

% For linearly spaced edges:
figure;
subplot(2,1,1);
plot(x,y,'.');
xedges = linspace(-1,1,64); yedges = linspace(-1,1,64);
histmat = hist2(x, y, xedges, yedges);
subplot(2,1,2);
% pcolor(xedges,yedges,histmat'); colorbar ; axis square tight ;

izero = histmat==0;
histmat(izero) = nan;
h=pcolor(xedges,yedges,histmat'); colorbar ; axis square tight ;
set(h,'edgecolor','none');
end

[Laszlo Balkay]

Comment to John: You should not transpose the histmat matrix:
  The pcolor(xedges,yedges,histmat) will work.

 

[Alwan]

Hi. How can I get the actual data points instead of its index. I actually want to add the magnitudes of the data that fall in each bin. Thank you