Help bining 2d data (or using hist2)

The following code does not seem to produce the right results. Any suggestions? The upper plot shows the points, while the lower plot shows the distribution of the points. The x and y ranges should be the same. Note, this example includes hist2 by Laszlo Balkay

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

function histmat  = hist2(x, y, xedges, yedges)
% 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.  
%
%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 ;

% University of Debrecen, PET Center/Laszlo Balkay 2006
% email: [e-mail]

if nargin ~= 4
    error ('The four input arguments are required!');
    return;
end
if any(size(x) ~= size(y)) 
    error ('The size of the two first input vectors should be same!');
    return;
end

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

%xbin, ybin zero for out of range values 
% (see the help of histc) force this event to the 
% first bins
xbin(find(xbin == 0)) = 1;
ybin(find(ybin == 0)) = 1;

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

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

%[xyuni, m, n] = unique(xy);
xyuni = unique(xy);
hstres = histc(xy,xyuni);
histmat = zeros(xnbin,ynbin);
histmat(xyuni-indexshift) = hstres;
histmat = histmat';

end

Answer 1

[John]

Never mind, the error is in the linspace range. Here's the correct code to show 2d hist.

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

end

% For linearly spaced edges:
figure;
subplot(2,1,1);
plot(x,y,'.');
% correct the linspace range
ax=axis;
xedges = linspace(ax(1),ax(2),64); yedges = linspace(ax(3),ax(4),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

function histmat  = hist2(x, y, xedges, yedges)
% 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.
%
%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 ;

% University of Debrecen, PET Center/Laszlo Balkay 2006
% email: [e-mail]

if nargin ~= 4
    error ('The four input arguments are required!');
    return;
end
if any(size(x) ~= size(y))
    error ('The size of the two first input vectors should be same!');
    return;
end

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

%xbin, ybin zero for out of range values
% (see the help of histc) force this event to the
% first bins
xbin((xbin == 0)) = 1;
ybin((ybin == 0)) = 1;

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

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

%[xyuni, m, n] = unique(xy);
xyuni = unique(xy);
hstres = histc(xy,xyuni);
histmat = zeros(xnbin,ynbin);
histmat(xyuni-indexshift) = hstres;
histmat = histmat';

end