Help bining 2d data (or using hist2)

3 views (last 30 days)
John
John on 1 Sep 2011
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: balkay@pet.dote.hu
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

Answers (1)

John
John on 1 Sep 2011
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: balkay@pet.dote.hu
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

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!